1//////////////////////////////////////////////////////////////////////
2// LibFile: geometry.scad
3// Perform calculations on lines, polygons, planes and circles, including
4// normals, intersections of objects, distance between objects, and tangent lines.
5// Throughout this library, lines can be treated as either unbounded lines, as rays with
6// a single endpoint or as segments, bounded by endpoints at both ends.
7// Includes:
8// include <BOSL2/std.scad>
9// FileGroup: Math
10// FileSummary: Geometrical calculations including intersections of lines, circles and planes, circle from 3 points
11// FileFootnotes: STD=Included in std.scad
12//////////////////////////////////////////////////////////////////////
13
14
15// Section: Lines, Rays, and Segments
16
17// Function: is_point_on_line()
18// Synopsis: Determine if a point is on a line, ray or segment.
19// Topics: Geometry, Points, Segments
20// See Also: is_collinear(), is_point_on_line(), point_line_distance(), line_from_points()
21// Usage:
22// pt = is_point_on_line(point, line, [bounded], [eps]);
23// Description:
24// Determine if the point is on the line segment, ray or segment defined by the two between two points.
25// Returns true if yes, and false if not. If bounded is set to true it specifies a segment, with
26// both lines bounded at the ends. Set bounded to `[true,false]` to get a ray. You can use
27// the shorthands RAY and SEGMENT to set bounded.
28// Arguments:
29// point = The point to test.
30// line = Array of two points defining the line, ray, or segment to test against.
31// bounded = boolean or list of two booleans defining endpoint conditions for the line. If false treat the line as an unbounded line. If true treat it as a segment. If [true,false] treat as a ray, based at the first endpoint. Default: false
32// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
33function is_point_on_line(point, line, bounded=false, eps=EPSILON) =
34 assert(is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
35 assert(is_vector(point), "Point must be a vector")
36 assert(_valid_line(line, len(point),eps),"Given line is not valid")
37 _is_point_on_line(point, line, bounded,eps);
38
39function _is_point_on_line(point, line, bounded=false, eps=EPSILON) =
40 let(
41 v1 = (line[1]-line[0]),
42 v0 = (point-line[0]),
43 t = v0*v1/(v1*v1),
44 bounded = force_list(bounded,2)
45 )
46 abs(cross(v0,v1))<=eps*norm(v1)
47 && (!bounded[0] || t>=-eps)
48 && (!bounded[1] || t<1+eps) ;
49
50
51///Internal - distance from point `d` to the line passing through the origin with unit direction n
52///_dist2line works for any dimension
53function _dist2line(d,n) = norm(d-(d * n) * n);
54
55
56///Internal
57function _valid_line(line,dim,eps=EPSILON) =
58 is_matrix(line,2,dim)
59 && norm(line[1]-line[0])>eps*max(norm(line[1]),norm(line[0]));
60
61//Internal
62function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps);
63
64
65/// Internal Function: _is_at_left()
66/// Usage:
67/// pt = point_left_of_line2d(point, line);
68/// Topics: Geometry, Points, Lines
69/// Description:
70/// Return true iff a 2d point is on or at left of the line defined by `line`.
71/// Arguments:
72/// pt = The 2d point to check position of.
73/// line = Array of two 2d points forming the line segment to test against.
74/// eps = Tolerance in the geometrical tests.
75function _is_at_left(pt,line,eps=EPSILON) = _tri_class([pt,line[0],line[1]],eps) <= 0;
76
77
78/// Internal Function: _degenerate_tri()
79/// Usage:
80/// degen = _degenerate_tri(triangle);
81/// Topics: Geometry, Triangles
82/// Description:
83/// Return true for a specific kind of degeneracy: any two triangle vertices are equal
84/// Arguments:
85/// tri = A list of three 2d points
86/// eps = Tolerance in the geometrical tests.
87function _degenerate_tri(tri,eps) =
88 max(norm(tri[0]-tri[1]), norm(tri[1]-tri[2]), norm(tri[2]-tri[0])) < eps ;
89
90
91/// Internal Function: _tri_class()
92/// Usage:
93/// class = _tri_class(triangle);
94/// Topics: Geometry, Triangles
95/// Description:
96/// Return 1 if the triangle `tri` is CW.
97/// Return 0 if the triangle `tri` has colinear vertices.
98/// Return -1 if the triangle `tri` is CCW.
99/// Arguments:
100/// tri = A list of the three 2d vertices of a triangle.
101/// eps = Tolerance in the geometrical tests.
102function _tri_class(tri, eps=EPSILON) =
103 let( crx = cross(tri[1]-tri[2],tri[0]-tri[2]) )
104 abs( crx ) <= eps*norm(tri[1]-tri[2])*norm(tri[0]-tri[2]) ? 0 : sign( crx );
105
106
107/// Internal Function: _pt_in_tri()
108/// Usage:
109/// class = _pt_in_tri(point, tri);
110/// Topics: Geometry, Points, Triangles
111/// Description:
112// For CW triangles `tri` :
113/// return 1 if point is inside the triangle interior.
114/// return =0 if point is on the triangle border.
115/// return -1 if point is outside the triangle.
116/// Arguments:
117/// point = The point to check position of.
118/// tri = A list of the three 2d vertices of a triangle.
119/// eps = Tolerance in the geometrical tests.
120function _pt_in_tri(point, tri, eps=EPSILON) =
121 min( _tri_class([tri[0],tri[1],point],eps),
122 _tri_class([tri[1],tri[2],point],eps),
123 _tri_class([tri[2],tri[0],point],eps) );
124
125
126/// Internal Function: _point_left_of_line2d()
127/// Usage:
128/// pt = point_left_of_line2d(point, line);
129/// Topics: Geometry, Points, Lines
130/// Description:
131/// Return >0 if point is left of the line defined by `line`.
132/// Return =0 if point is on the line.
133/// Return <0 if point is right of the line.
134/// Arguments:
135/// point = The point to check position of.
136/// line = Array of two points forming the line segment to test against.
137function _point_left_of_line2d(point, line, eps=EPSILON) =
138 assert( is_vector(point,2) && is_vector(line*point, 2), "Improper input." )
139// cross(line[0]-point, line[1]-line[0]);
140 _tri_class([point,line[1],line[0]],eps);
141
142
143// Function: is_collinear()
144// Synopsis: Determine if points are collinear.
145// Topics: Geometry, Points, Collinearity
146// See Also: is_collinear(), is_point_on_line(), point_line_distance(), line_from_points()
147// Usage:
148// bool = is_collinear(a, [b, c], [eps]);
149// Description:
150// Returns true if the points `a`, `b` and `c` are co-linear or if the list of points `a` is collinear.
151// Arguments:
152// a = First point or list of points.
153// b = Second point or undef; it should be undef if `c` is undef
154// c = Third point or undef.
155// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
156function is_collinear(a, b, c, eps=EPSILON) =
157 assert( is_path([a,b,c],dim=undef)
158 || ( is_undef(b) && is_undef(c) && is_path(a,dim=undef) ),
159 "Input should be 3 points or a list of points with same dimension.")
160 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
161 let( points = is_def(c) ? [a,b,c]: a )
162 len(points)<3 ? true :
163 _noncollinear_triple(points,error=false,eps=eps) == [];
164
165
166// Function: point_line_distance()
167// Synopsis: Find shortest distance from point to a line, segment or ray.
168// Topics: Geometry, Points, Lines, Distance
169// See Also: is_collinear(), is_point_on_line(), point_line_distance(), line_from_points()
170// Usage:
171// dist = point_line_distance(pt, line, [bounded]);
172// Description:
173// Finds the shortest distance from the point `pt` to the specified line, segment or ray.
174// The bounded parameter specifies the whether the endpoints give a ray or segment.
175// By default assumes an unbounded line.
176// Arguments:
177// pt = A point to find the distance of from the line.
178// line = A list of two points defining a line.
179// bounded = a boolean or list of two booleans specifiying whether each end is bounded. Default: false
180// Example:
181// dist1 = point_line_distance([3,8], [[-10,0], [10,0]]); // Returns: 8
182// dist2 = point_line_distance([3,8], [[-10,0], [10,0]],SEGMENT); // Returns: 8
183// dist3 = point_line_distance([14,3], [[-10,0], [10,0]],SEGMENT); // Returns: 5
184function point_line_distance(pt, line, bounded=false) =
185 assert(is_bool(bounded) || is_bool_list(bounded,2), "\"bounded\" is invalid")
186 assert( _valid_line(line) && is_vector(pt,len(line[0])),
187 "Invalid line, invalid point or incompatible dimensions." )
188 bounded == LINE ? _dist2line(pt-line[0],unit(line[1]-line[0]))
189 : norm(pt-line_closest_point(line,pt,bounded));
190
191
192// Function: segment_distance()
193// Synopsis: Find smallest distance between two line semgnets.
194// Topics: Geometry, Segments, Distance
195// See Also: convex_collision(), convex_distance()
196// Usage:
197// dist = segment_distance(seg1, seg2, [eps]);
198// Description:
199// Returns the smallest distance of the points on two given line segments.
200// Arguments:
201// seg1 = The list of two points representing the first line segment to check the distance of.
202// seg2 = The list of two points representing the second line segment to check the distance of.
203// eps = tolerance for point comparisons
204// Example:
205// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5
206// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0
207function segment_distance(seg1, seg2,eps=EPSILON) =
208 assert( is_matrix(concat(seg1,seg2),4), "Inputs should be two valid segments." )
209 convex_distance(seg1,seg2,eps);
210
211
212// Function: line_normal()
213// Synopsis: Return normal vector to given line.
214// Topics: Geometry, Lines
215// See Also: line_intersection(), line_from_points()
216// Usage:
217// vec = line_normal([P1,P2])
218// vec = line_normal(p1,p2)
219// Description:
220// Returns the 2D normal vector to the given 2D line. This is otherwise known as the perpendicular vector counter-clockwise to the given ray.
221// Arguments:
222// p1 = First point on 2D line.
223// p2 = Second point on 2D line.
224// Example(2D):
225// p1 = [10,10];
226// p2 = [50,30];
227// n = line_normal(p1,p2);
228// stroke([p1,p2], endcap2="arrow2");
229// color("green") stroke([p1,p1+10*n], endcap2="arrow2");
230// color("blue") move_copies([p1,p2]) circle(d=2, $fn=12);
231function line_normal(p1,p2) =
232 is_undef(p2)
233 ? assert( len(p1)==2 && !is_undef(p1[1]) , "Invalid input." )
234 line_normal(p1[0],p1[1])
235 : assert( _valid_line([p1,p2],dim=2), "Invalid line." )
236 unit([p1.y-p2.y,p2.x-p1.x]);
237
238
239// 2D Line intersection from two segments.
240// This function returns [p,t,u] where p is the intersection point of
241// the lines defined by the two segments, t is the proportional distance
242// of the intersection point along s1, and u is the proportional distance
243// of the intersection point along s2. The proportional values run over
244// the range of 0 to 1 for each segment, so if it is in this range, then
245// the intersection lies on the segment. Otherwise it lies somewhere on
246// the extension of the segment. If lines are parallel or coincident then
247// it returns undef.
248
249function _general_line_intersection(s1,s2,eps=EPSILON) =
250 let(
251 denominator = cross(s1[0]-s1[1],s2[0]-s2[1])
252 )
253 approx(denominator,0,eps=eps) ? undef :
254 let(
255 t = cross(s1[0]-s2[0],s2[0]-s2[1]) / denominator,
256 u = cross(s1[0]-s2[0],s1[0]-s1[1]) / denominator
257 )
258 [s1[0]+t*(s1[1]-s1[0]), t, u];
259
260
261// Function: line_intersection()
262// Synopsis: Compute intersection of two lines, segments or rays.
263// Topics: Geometry, Lines
264// See Also: line_normal(), line_from_points()
265// Usage:
266// pt = line_intersection(line1, line2, [bounded1], [bounded2], [bounded=], [eps=]);
267// Description:
268// Returns the intersection point of any two 2D lines, segments or rays. Returns undef
269// if they do not intersect. You specify a line by giving two distinct points on the
270// line. You specify rays or segments by giving a pair of points and indicating
271// bounded[0]=true to bound the line at the first point, creating rays based at l1[0] and l2[0],
272// or bounded[1]=true to bound the line at the second point, creating the reverse rays bounded
273// at l1[1] and l2[1]. If bounded=[true, true] then you have segments defined by their two
274// endpoints. By using bounded1 and bounded2 you can mix segments, rays, and lines as needed.
275// You can set the bounds parameters to true as a shorthand for [true,true] to sepcify segments.
276// Arguments:
277// line1 = List of two points in 2D defining the first line, segment or ray
278// line2 = List of two points in 2D defining the second line, segment or ray
279// bounded1 = boolean or list of two booleans defining which ends are bounded for line1. Default: [false,false]
280// bounded2 = boolean or list of two booleans defining which ends are bounded for line2. Default: [false,false]
281// ---
282// bounded = boolean or list of two booleans defining which ends are bounded for both lines. The bounded1 and bounded2 parameters override this if both are given.
283// eps = tolerance for geometric comparisons. Default: `EPSILON` (1e-9)
284// Example(2D): The segments do not intersect but the lines do in this example.
285// line1 = 10*[[9, 4], [5, 7]];
286// line2 = 10*[[2, 3], [6, 5]];
287// stroke(line1, endcaps="arrow2");
288// stroke(line2, endcaps="arrow2");
289// isect = line_intersection(line1, line2);
290// color("red") translate(isect) circle(r=1,$fn=12);
291// Example(2D): Specifying a ray and segment using the shorthand variables.
292// line1 = 10*[[0, 2], [4, 7]];
293// line2 = 10*[[10, 4], [3, 4]];
294// stroke(line1);
295// stroke(line2, endcap2="arrow2");
296// isect = line_intersection(line1, line2, SEGMENT, RAY);
297// color("red") translate(isect) circle(r=1,$fn=12);
298// Example(2D): Here we use the same example as above, but specify two segments using the bounded argument.
299// line1 = 10*[[0, 2], [4, 7]];
300// line2 = 10*[[10, 4], [3, 4]];
301// stroke(line1);
302// stroke(line2);
303// isect = line_intersection(line1, line2, bounded=true); // Returns undef
304function line_intersection(line1, line2, bounded1, bounded2, bounded, eps=EPSILON) =
305 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
306 assert( _valid_line(line1,dim=2,eps=eps), "First line invalid")
307 assert( _valid_line(line2,dim=2,eps=eps), "Second line invalid")
308 assert( is_undef(bounded) || is_bool(bounded) || is_bool_list(bounded,2), "Invalid value for \"bounded\"")
309 assert( is_undef(bounded1) || is_bool(bounded1) || is_bool_list(bounded1,2), "Invalid value for \"bounded1\"")
310 assert( is_undef(bounded2) || is_bool(bounded2) || is_bool_list(bounded2,2), "Invalid value for \"bounded2\"")
311 let(isect = _general_line_intersection(line1,line2,eps=eps))
312 is_undef(isect) ? undef :
313 let(
314 bounded1 = force_list(first_defined([bounded1,bounded,false]),2),
315 bounded2 = force_list(first_defined([bounded2,bounded,false]),2),
316 good = (!bounded1[0] || isect[1]>=0-eps)
317 && (!bounded1[1] || isect[1]<=1+eps)
318 && (!bounded2[0] || isect[2]>=0-eps)
319 && (!bounded2[1] || isect[2]<=1+eps)
320 )
321 good ? isect[0] : undef;
322
323
324// Function: line_closest_point()
325// Synopsis: Find point on given line, segment or ray that is closest to a given point.
326// Topics: Geometry, Lines, Distance
327// See Also: line_normal(), point_line_distance()
328// Usage:
329// pt = line_closest_point(line, pt, [bounded]);
330// Description:
331// Returns the point on the given line, segment or ray that is closest to the given point `pt`.
332// The inputs `line` and `pt` args should be of the same dimension. The parameter bounded indicates
333// whether the points of `line` should be treated as endpoints.
334// Arguments:
335// line = A list of two points that are on the unbounded line.
336// pt = The point to find the closest point on the line to.
337// bounded = boolean or list of two booleans indicating that the line is bounded at that end. Default: [false,false]
338// Example(2D):
339// line = [[-30,0],[30,30]];
340// pt = [-32,-10];
341// p2 = line_closest_point(line,pt);
342// stroke(line, endcaps="arrow2");
343// color("blue") translate(pt) circle(r=1,$fn=12);
344// color("red") translate(p2) circle(r=1,$fn=12);
345// Example(2D): If the line is bounded on the left you get the endpoint instead
346// line = [[-30,0],[30,30]];
347// pt = [-32,-10];
348// p2 = line_closest_point(line,pt,bounded=[true,false]);
349// stroke(line, endcap2="arrow2");
350// color("blue") translate(pt) circle(r=1,$fn=12);
351// color("red") translate(p2) circle(r=1,$fn=12);
352// Example(2D): In this case it doesn't matter how bounded is set. Using SEGMENT is the most restrictive option.
353// line = [[-30,0],[30,30]];
354// pt = [-5,0];
355// p2 = line_closest_point(line,pt,SEGMENT);
356// stroke(line);
357// color("blue") translate(pt) circle(r=1,$fn=12);
358// color("red") translate(p2) circle(r=1,$fn=12);
359// Example(2D): The result here is the same for a line or a ray.
360// line = [[-30,0],[30,30]];
361// pt = [40,25];
362// p2 = line_closest_point(line,pt,RAY);
363// stroke(line, endcap2="arrow2");
364// color("blue") translate(pt) circle(r=1,$fn=12);
365// color("red") translate(p2) circle(r=1,$fn=12);
366// Example(2D): But with a segment we get a different result
367// line = [[-30,0],[30,30]];
368// pt = [40,25];
369// p2 = line_closest_point(line,pt,SEGMENT);
370// stroke(line);
371// color("blue") translate(pt) circle(r=1,$fn=12);
372// color("red") translate(p2) circle(r=1,$fn=12);
373// Example(2D): The shorthand RAY uses the first point as the base of the ray. But you can specify a reversed ray directly, and in this case the result is the same as the result above for the segment.
374// line = [[-30,0],[30,30]];
375// pt = [40,25];
376// p2 = line_closest_point(line,pt,[false,true]);
377// stroke(line,endcap1="arrow2");
378// color("blue") translate(pt) circle(r=1,$fn=12);
379// color("red") translate(p2) circle(r=1,$fn=12);
380// Example(FlatSpin,VPD=200,VPT=[0,0,15]): A 3D example
381// line = [[-30,-15,0],[30,15,30]];
382// pt = [5,5,5];
383// p2 = line_closest_point(line,pt);
384// stroke(line, endcaps="arrow2");
385// color("blue") translate(pt) sphere(r=1,$fn=12);
386// color("red") translate(p2) sphere(r=1,$fn=12);
387function line_closest_point(line, pt, bounded=false) =
388 assert(_valid_line(line), "Invalid line")
389 assert(is_vector(pt, len(line[0])), "Invalid point or incompatible dimensions.")
390 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid value for \"bounded\"")
391 let(
392 bounded = force_list(bounded,2)
393 )
394 bounded==[false,false] ?
395 let( n = unit( line[0]- line[1]) )
396 line[1] + ((pt- line[1]) * n) * n
397 : bounded == [true,true] ?
398 pt + _closest_s1([line[0]-pt, line[1]-pt])[0]
399 :
400 let(
401 ray = bounded==[true,false] ? line : reverse(line),
402 seglen = norm(ray[1]-ray[0]),
403 segvec = (ray[1]-ray[0])/seglen,
404 projection = (pt-ray[0]) * segvec
405 )
406 projection<=0 ? ray[0] :
407 ray[0] + projection*segvec;
408
409
410// Function: line_from_points()
411// Synopsis: Given a list of collinear points, return the line they define.
412// Topics: Geometry, Lines, Points
413// Usage:
414// line = line_from_points(points, [fast], [eps]);
415// Description:
416// Given a list of 2 or more collinear points, returns two points defining a line containing them.
417// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned.
418// if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned.
419// Arguments:
420// points = The list of points to find the line through.
421// fast = If true, don't verify that all points are collinear. Default: false
422// eps = How much variance is allowed in testing each point against the line. Default: `EPSILON` (1e-9)
423function line_from_points(points, fast=false, eps=EPSILON) =
424 assert( is_path(points), "Invalid point list." )
425 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
426 let( pb = furthest_point(points[0],points) )
427 norm(points[pb]-points[0])<eps*max(norm(points[pb]),norm(points[0])) ? undef :
428 fast || is_collinear(points)
429 ? [points[pb], points[0]]
430 : undef;
431
432
433
434// Section: Planes
435
436
437// Function: is_coplanar()
438// Synopsis: Check if 3d points are coplanar and not collinear.
439// Topics: Geometry, Coplanarity
440// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon()
441// Usage:
442// bool = is_coplanar(points,[eps]);
443// Description:
444// Returns true if the given 3D points are non-collinear and are on a plane.
445// Arguments:
446// points = The points to test.
447// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
448function is_coplanar(points, eps=EPSILON) =
449 assert( is_path(points,dim=3) , "Input should be a list of 3D points." )
450 assert( is_finite(eps) && eps>=0, "The tolerance should be a non-negative value." )
451 len(points)<=2 ? false
452 : let( ip = _noncollinear_triple(points,error=false,eps=eps) )
453 ip == [] ? false :
454 let( plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]) )
455 _pointlist_greatest_distance(points,plane) < eps;
456
457
458
459// Function: plane3pt()
460// Synopsis: Return a plane from 3 points.
461// Topics: Geometry, Planes
462// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon()
463// Usage:
464// plane = plane3pt(p1, p2, p3);
465// plane = plane3pt([p1, p2, p3]);
466// Description:
467// Generates the normalized cartesian equation of a plane from three 3d points.
468// Returns [A,B,C,D] where Ax + By + Cz = D is the equation of a plane.
469// Returns undef, if the points are collinear.
470// Arguments:
471// p1 = The first point on the plane.
472// p2 = The second point on the plane.
473// p3 = The third point on the plane.
474function plane3pt(p1, p2, p3) =
475 is_undef(p2) && is_undef(p3) && is_path(p1,dim=3) ? plane3pt(p1[0],p1[1],p1[2])
476 : assert( is_path([p1,p2,p3],dim=3) && len(p1)==3,
477 "Invalid points or incompatible dimensions." )
478 let(
479 crx = cross(p3-p1, p2-p1),
480 nrm = norm(crx)
481 ) approx(nrm,0) ? undef :
482 concat(crx, crx*p1)/nrm;
483
484
485// Function: plane3pt_indexed()
486// Synopsis: Given list of 3d points and 3 indices, return the plane they define.
487// Topics: Geometry, Planes
488// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon()
489// Usage:
490// plane = plane3pt_indexed(points, i1, i2, i3);
491// Description:
492// Given a list of 3d points, and the indices of three of those points,
493// generates the normalized cartesian equation of a plane that those points all
494// lie on. If the points are not collinear, returns [A,B,C,D] where Ax+By+Cz=D is the equation of a plane.
495// If they are collinear, returns [].
496// Arguments:
497// points = A list of points.
498// i1 = The index into `points` of the first point on the plane.
499// i2 = The index into `points` of the second point on the plane.
500// i3 = The index into `points` of the third point on the plane.
501function plane3pt_indexed(points, i1, i2, i3) =
502 is_undef(i3) && is_undef(i2) && is_vector(i1) ? plane3pt_indexed(points, i1[0], i1[1], i1[2])
503 :
504 assert( is_vector([i1,i2,i3]) && min(i1,i2,i3)>=0 && is_list(points) && max(i1,i2,i3)<len(points),
505 "Invalid or out of range indices." )
506 assert( is_path([points[i1], points[i2], points[i3]],dim=3),
507 "Improper points or improper dimensions." )
508 let(
509 p1 = points[i1],
510 p2 = points[i2],
511 p3 = points[i3]
512 ) plane3pt(p1,p2,p3);
513
514
515// Function: plane_from_normal()
516// Synopsis: Return plane defined by normal vector and a point.
517// Topics: Geometry, Planes
518// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon()
519// Usage:
520// plane = plane_from_normal(normal, [pt])
521// Description:
522// Returns a plane defined by a normal vector and a point. If you omit `pt` you will get a plane
523// passing through the origin.
524// Arguments:
525// normal = Normal vector to the plane to find.
526// pt = Point 3D on the plane to find.
527// Example:
528// plane_from_normal([0,0,1], [2,2,2]); // Returns the xy plane passing through the point (2,2,2)
529function plane_from_normal(normal, pt=[0,0,0]) =
530 assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),0),
531 "Inputs `normal` and `pt` should be 3d vectors/points and `normal` cannot be zero." )
532 concat(normal, normal*pt) / norm(normal);
533
534
535// Eigenvalues for a 3x3 symmetrical matrix in decreasing order
536// Based on: https://en.wikipedia.org/wiki/Eigenvalue_algorithm
537function _eigenvals_symm_3(M) =
538 let( p1 = pow(M[0][1],2) + pow(M[0][2],2) + pow(M[1][2],2) )
539 (p1<EPSILON)
540 ? -sort(-[ M[0][0], M[1][1], M[2][2] ]) // diagonal matrix: eigenvals in decreasing order
541 : let( q = (M[0][0]+M[1][1]+M[2][2])/3,
542 B = (M - q*ident(3)),
543 dB = [B[0][0], B[1][1], B[2][2]],
544 p2 = dB*dB + 2*p1,
545 p = sqrt(p2/6),
546 r = det3(B/p)/2,
547 ph = acos(constrain(r,-1,1))/3,
548 e1 = q + 2*p*cos(ph),
549 e3 = q + 2*p*cos(ph+120),
550 e2 = 3*q - e1 - e3 )
551 [ e1, e2, e3 ];
552
553
554// the i-th normalized eigenvector of a 3x3 symmetrical matrix M from its eigenvalues
555// using Cayley–Hamilton theorem according to:
556// https://en.wikipedia.org/wiki/Eigenvalue_algorithm
557function _eigenvec_symm_3(M,evals,i=0) =
558 let(
559 I = ident(3),
560 A = (M - evals[(i+1)%3]*I) * (M - evals[(i+2)%3]*I) ,
561 k = max_index( [for(i=[0:2]) norm(A[i]) ])
562 )
563 norm(A[k])<EPSILON ? I[k] : A[k]/norm(A[k]);
564
565
566// finds the eigenvector corresponding to the smallest eigenvalue of the covariance matrix of a pointlist
567// returns the mean of the points, the eigenvector and the greatest eigenvalue
568function _covariance_evec_eval(points) =
569 let( pm = sum(points)/len(points), // mean point
570 Y = [ for(i=[0:len(points)-1]) points[i] - pm ],
571 M = transpose(Y)*Y , // covariance matrix
572 evals = _eigenvals_symm_3(M), // eigenvalues in decreasing order
573 evec = _eigenvec_symm_3(M,evals,i=2) )
574 [pm, evec, evals[0] ];
575
576
577// Function: plane_from_points()
578// Synopsis: Return plane defined by a set of coplanar 3d points, with arbitrary normal direction.
579// Topics: Geometry, Planes, Points
580// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon()
581// Usage:
582// plane = plane_from_points(points, [fast], [eps]);
583// Description:
584// Given a list of 3 or more coplanar 3D points, returns the coefficients of the normalized cartesian equation of a plane,
585// that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane and norm([A,B,C])=1.
586// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned.
587// If `fast` is true, the polygon coplanarity check is skipped and a best fitting plane is returned.
588// The direction of the plane's normal is arbitrary and is not determined by the point order, unlike {{plane_from_polygon()}}.
589// This function is faster than {{plane_from_polygon()}}.
590// Arguments:
591// points = The list of points to find the plane of.
592// fast = If true, don't verify the point coplanarity. Default: false
593// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
594// Example(3D):
595// points = rot(45, v=[-0.3,1,0], p=path3d(random_points(25,2,scale=55,seed=47), 70));
596// plane = plane_from_points(points);
597// #move_copies(points)sphere(d=3);
598// cp = mean(points);
599// move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(50);
600function plane_from_points(points, fast=false, eps=EPSILON) =
601 assert( is_path(points,dim=3), "Improper 3d point list." )
602 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
603 len(points) == 3
604 ? plane3pt(points[0],points[1],points[2])
605 : let(
606 covmix = _covariance_evec_eval(points),
607 pm = covmix[0],
608 evec = covmix[1],
609 eval0 = covmix[2],
610 plane = [ each evec, pm*evec]
611 )
612 !fast && _pointlist_greatest_distance(points,plane)>eps*eval0 ? undef :
613 plane ;
614
615
616// Function: plane_from_polygon()
617// Synopsis: Given a 3d planar polygon, returns directed plane.
618// Topics: Geometry, Planes, Polygons
619// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon()
620// Usage:
621// plane = plane_from_polygon(points, [fast], [eps]);
622// Description:
623// Given a 3D planar polygon, returns the normalized cartesian equation of its plane.
624// Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1.
625// If not all the points in the polygon are coplanar, then [] is returned.
626// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned.
627// if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned.
628// The normal direction is determined by the order of the points and the right hand rule. This is slower than {{plane_from_points()}},
629// which returns an arbitrary normal.
630// Arguments:
631// poly = The planar 3D polygon to find the plane of.
632// fast = If true, doesn't verify that all points in the polygon are coplanar. Default: false
633// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
634// Example(3D):
635// xyzpath = rot(45, v=[0,1,0], p=path3d(star(n=5,step=2,d=100), 70));
636// plane = plane_from_polygon(xyzpath);
637// #stroke(xyzpath,closed=true,width=3);
638// cp = centroid(xyzpath);
639// move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(45);
640function plane_from_polygon(poly, fast=false, eps=EPSILON) =
641 assert( is_path(poly,dim=3), "Invalid polygon." )
642 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
643 let(
644 poly_normal = polygon_normal(poly)
645 )
646 is_undef(poly_normal) ? undef :
647 let(
648 plane = plane_from_normal(poly_normal, poly[0])
649 )
650 fast? plane: are_points_on_plane(poly, plane, eps=eps)? plane: undef;
651
652
653// Function: plane_normal()
654// Synopsis: Returns the normal vector to a plane.
655// Topics: Geometry, Planes
656// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon(), plane_normal(), plane_offset()
657// Usage:
658// vec = plane_normal(plane);
659// Description:
660// Returns the unit length normal vector for the given plane.
661// Arguments:
662// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
663function plane_normal(plane) =
664 assert( _valid_plane(plane), "Invalid input plane." )
665 unit([plane.x, plane.y, plane.z]);
666
667
668// Function: plane_offset()
669// Synopsis: Returns the signed offset of the plane from the origin.
670// Topics: Geometry, Planes
671// See Also: plane3pt(), plane3pt_indexed(), plane_from_normal(), plane_from_points(), plane_from_polygon(), plane_normal(), plane_offset()
672// Usage:
673// d = plane_offset(plane);
674// Description:
675// Returns coeficient D of the normalized plane equation `Ax+By+Cz=D`, or the scalar offset of the plane from the origin.
676// This value may be negative.
677// The absolute value of this coefficient is the distance of the plane from the origin.
678// Arguments:
679// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
680function plane_offset(plane) =
681 assert( _valid_plane(plane), "Invalid input plane." )
682 plane[3]/norm([plane.x, plane.y, plane.z]);
683
684
685
686// Returns [POINT, U] if line intersects plane at one point, where U is zero at line[0] and 1 at line[1]
687// Returns [LINE, undef] if the line is on the plane.
688// Returns undef if line is parallel to, but not on the given plane.
689function _general_plane_line_intersection(plane, line, eps=EPSILON) =
690 let(
691 a = plane*[each line[0],-1], // evaluation of the plane expression at line[0]
692 b = plane*[each(line[1]-line[0]),0] // difference between the plane expression evaluation at line[1] and at line[0]
693 )
694 approx(b,0,eps) // is (line[1]-line[0]) "parallel" to the plane ?
695 ? approx(a,0,eps) // is line[0] on the plane ?
696 ? [line,undef] // line is on the plane
697 : undef // line is parallel but not on the plane
698 : [ line[0]-a/b*(line[1]-line[0]), -a/b ];
699
700
701/// Internal Function: normalize_plane()
702/// Usage:
703/// nplane = normalize_plane(plane);
704/// Topics: Geometry, Planes
705/// Description:
706/// Returns a new representation [A,B,C,D] of `plane` where norm([A,B,C]) is equal to one.
707function _normalize_plane(plane) =
708 assert( _valid_plane(plane), str("Invalid plane. ",plane ) )
709 plane/norm(point3d(plane));
710
711
712// Function: plane_line_intersection()
713// Synopsis: Returns the intersection of a plane and 3d line, segment or ray.
714// Topics: Geometry, Planes, Lines, Intersection
715// See Also: plane3pt(), plane_from_normal(), plane_from_points(), plane_from_polygon(), line_intersection()
716// Usage:
717// pt = plane_line_intersection(plane, line, [bounded], [eps]);
718// Description:
719// Takes a line, and a plane [A,B,C,D] where the equation of that plane is `Ax+By+Cz=D`.
720// If `line` intersects `plane` at one point, then that intersection point is returned.
721// If `line` lies on `plane`, then the original given `line` is returned.
722// If `line` is parallel to, but not on `plane`, then undef is returned.
723// Arguments:
724// plane = The [A,B,C,D] values for the equation of the plane.
725// line = A list of two distinct 3D points that are on the line.
726// bounded = If false, the line is considered unbounded. If true, it is treated as a bounded line segment. If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray. Default: false (unbounded)
727// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
728function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
729 assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
730 assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.")
731 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.")
732 let(
733 bounded = is_list(bounded)? bounded : [bounded, bounded],
734 res = _general_plane_line_intersection(plane, line, eps=eps)
735 ) is_undef(res) ? undef :
736 is_undef(res[1]) ? res[0] :
737 bounded[0] && res[1]<0 ? undef :
738 bounded[1] && res[1]>1 ? undef :
739 res[0];
740
741
742
743// Function: plane_intersection()
744// Synopsis: Returns the intersection of two or three planes.
745// Topics: Geometry, Planes, Intersection
746// See Also: plane3pt(), plane_from_normal(), plane_from_points(), plane_from_polygon(), line_intersection()
747// Usage:
748// line = plane_intersection(plane1, plane2)
749// pt = plane_intersection(plane1, plane2, plane3)
750// Description:
751// Compute the point which is the intersection of the three planes, or the line intersection of two planes.
752// If you give three planes the intersection is returned as a point. If you give two planes the intersection
753// is returned as a list of two points on the line of intersection. If any two input planes are parallel
754// or coincident then returns undef.
755// Arguments:
756// plane1 = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`.
757// plane2 = The [A,B,C,D] coefficients for the second plane equation `Ax+By+Cz=D`.
758// plane3 = The [A,B,C,D] coefficients for the third plane equation `Ax+By+Cz=D`.
759function plane_intersection(plane1,plane2,plane3) =
760 assert( _valid_plane(plane1) && _valid_plane(plane2) && (is_undef(plane3) ||_valid_plane(plane3)),
761 "The input must be 2 or 3 planes." )
762 is_def(plane3)
763 ? let(
764 matrix = [for(p=[plane1,plane2,plane3]) point3d(p)],
765 rhs = [for(p=[plane1,plane2,plane3]) p[3]]
766 )
767 linear_solve(matrix,rhs)
768 : let( normal = cross(plane_normal(plane1), plane_normal(plane2)) )
769 approx(norm(normal),0) ? undef :
770 let(
771 matrix = [for(p=[plane1,plane2]) point3d(p)],
772 rhs = [plane1[3], plane2[3]],
773 point = linear_solve(matrix,rhs)
774 )
775 point==[]? undef:
776 [point, point+normal];
777
778
779
780// Function: plane_line_angle()
781// Synopsis: Returns the angle between a plane and a 3d line.
782// Topics: Geometry, Planes, Lines, Angle
783// See Also: plane3pt(), plane_from_normal(), plane_from_points(), plane_from_polygon(), plane_intersection(), line_intersection(), vector_angle()
784// Usage:
785// angle = plane_line_angle(plane,line);
786// Description:
787// Compute the angle between a plane [A, B, C, D] and a 3d line, specified as a pair of 3d points [p1,p2].
788// The resulting angle is signed, with the sign positive if the vector p2-p1 lies above the plane, on
789// the same side of the plane as the plane's normal vector.
790function plane_line_angle(plane, line) =
791 assert( _valid_plane(plane), "Invalid plane." )
792 assert( _valid_line(line,dim=3), "Invalid 3d line." )
793 let(
794 linedir = unit(line[1]-line[0]),
795 normal = plane_normal(plane),
796 sin_angle = linedir*normal,
797 cos_angle = norm(cross(linedir,normal))
798 ) atan2(sin_angle,cos_angle);
799
800
801
802// Function: plane_closest_point()
803// Synopsis: Returns the orthogonal projection of points onto a plane.
804// Topics: Geometry, Planes, Projection
805// See Also: plane3pt(), line_closest_point(), point_plane_distance()
806// Usage:
807// pts = plane_closest_point(plane, points);
808// Description:
809// Given a plane definition `[A,B,C,D]`, where `Ax+By+Cz=D`, and a list of 2d or
810// 3d points, return the closest 3D orthogonal projection of the points on the plane.
811// In other words, for every point given, returns the closest point to it on the plane.
812// If points is a single point then returns a single point result.
813// Arguments:
814// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
815// points = List of points to project
816// Example(FlatSpin,VPD=500,VPT=[2,20,10]):
817// points = move([10,20,30], p=yrot(25, p=path3d(circle(d=100, $fn=36))));
818// plane = plane_from_normal([1,0,1]);
819// proj = plane_closest_point(plane,points);
820// color("red") move_copies(points) sphere(d=4,$fn=12);
821// color("blue") move_copies(proj) sphere(d=4,$fn=12);
822// move(centroid(proj)) {
823// rot(from=UP,to=plane_normal(plane)) {
824// anchor_arrow(50);
825// %cube([120,150,0.1],center=true);
826// }
827// }
828function plane_closest_point(plane, points) =
829 is_vector(points,3) ? plane_closest_point(plane,[points])[0] :
830 assert( _valid_plane(plane), "Invalid plane." )
831 assert( is_matrix(points,undef,3), "Must supply 3D points.")
832 let(
833 plane = _normalize_plane(plane),
834 n = point3d(plane)
835 )
836 [for(pi=points) pi - (pi*n - plane[3])*n];
837
838
839// Function: point_plane_distance()
840// Synopsis: Determine distance between a point and plane.
841// Topics: Geometry, Planes, Distance
842// See Also: plane3pt(), line_closest_point(), plane_closest_point()
843// Usage:
844// dist = point_plane_distance(plane, point)
845// Description:
846// Given a plane as [A,B,C,D] where the cartesian equation for that plane
847// is Ax+By+Cz=D, determines how far from that plane the given point is.
848// The returned distance will be positive if the point is above the
849// plane, meaning on the side where the plane normal points.
850// If the point is below the plane, then the distance returned
851// will be negative. The normal of the plane is [A,B,C].
852// Arguments:
853// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
854// point = The distance evaluation point.
855function point_plane_distance(plane, point) =
856 assert( _valid_plane(plane), "Invalid input plane." )
857 assert( is_vector(point,3), "The point should be a 3D point." )
858 let( plane = _normalize_plane(plane) )
859 point3d(plane)* point - plane[3];
860
861
862
863// the maximum distance from points to the plane
864function _pointlist_greatest_distance(points,plane) =
865 let(
866 normal = [plane[0],plane[1],plane[2]],
867 pt_nrm = points*normal
868 )
869 max( max(pt_nrm) - plane[3], -min(pt_nrm) + plane[3]) / norm(normal);
870
871
872// Function: are_points_on_plane()
873// Synopsis: Determine if all of the listed points are on a plane.
874// Topics: Geometry, Planes, Points
875// See Also: plane3pt(), line_closest_point(), plane_closest_point(), is_coplanar()
876// Usage:
877// bool = are_points_on_plane(points, plane, [eps]);
878// Description:
879// Returns true if the given 3D points are on the given plane.
880// Arguments:
881// plane = The plane to test the points on.
882// points = The list of 3D points to test.
883// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
884function are_points_on_plane(points, plane, eps=EPSILON) =
885 assert( _valid_plane(plane), "Invalid plane." )
886 assert( is_matrix(points,undef,3) && len(points)>0, "Invalid pointlist." ) // using is_matrix it accepts len(points)==1
887 assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
888 _pointlist_greatest_distance(points,plane) < eps;
889
890
891/// Internal Function: is_point_above_plane()
892/// Usage:
893/// bool = _is_point_above_plane(plane, point);
894/// Topics: Geometry, Planes
895/// Description:
896/// Given a plane as [A,B,C,D] where the cartesian equation for that plane
897/// is Ax+By+Cz=D, determines if the given 3D point is on the side of that
898/// plane that the normal points towards. The normal of the plane is the
899/// same as [A,B,C].
900/// Arguments:
901/// plane = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`.
902/// point = The 3D point to test.
903function _is_point_above_plane(plane, point) =
904 point_plane_distance(plane, point) > EPSILON;
905
906
907
908// Section: Circle Calculations
909
910// Function: circle_line_intersection()
911// Synopsis: Find the intersection points between a 2d circle and a line, ray or segment.
912// Topics: Geometry, Circles, Lines, Intersection
913// See Also: circle_line_intersection(), circle_circle_intersection(), circle_2tangents(), circle_3points(), circle_point_tangents(), circle_circle_tangents()
914// Usage:
915// pts = circle_line_intersection(r|d=, cp, line, [bounded], [eps=]);
916// Description:
917// Find intersection points between a 2D circle and a line, ray or segment specified by two points.
918// By default the line is unbounded. Returns the list of zero or more intersection points.
919// Arguments:
920// r = Radius of circle
921// cp = Center of circle
922// line = Two points defining the line
923// bounded = False for unbounded line, true for a segment, or a vector [false,true] or [true,false] to specify a ray with the first or second end unbounded. Default: false
924// ---
925// d = Diameter of circle
926// eps = Epsilon used for identifying the case with one solution. Default: `1e-9`
927// Example(2D): Standard intersection returns two points.
928// line = [[-15,2], [15,7]];
929// cp = [1,2]; r = 10;
930// translate(cp) circle(r=r);
931// color("black") stroke(line, endcaps="arrow2", width=0.5);
932// isects = circle_line_intersection(r=r, cp=cp, line=line);
933// color("red") move_copies(isects) circle(d=1);
934// Example(2D): Tangent intersection returns one point.
935// line = [[-10,12], [10,12]];
936// cp = [1,2]; r = 10;
937// translate(cp) circle(r=r);
938// color("black") stroke(line, endcaps="arrow2", width=0.5);
939// isects = circle_line_intersection(r=r, cp=cp, line=line);
940// color("#f44") move_copies(isects) circle(d=1);
941// Example(2D): A bounded ray might only intersect in one direction.
942// line = [[-5,2], [5,7]];
943// extended = [line[0], line[0]+22*unit(line[1]-line[0])];
944// cp = [1,2]; r = 10;
945// translate(cp) circle(r=r);
946// color("gray") dashed_stroke(extended, width=0.2);
947// color("black") stroke(line, endcap2="arrow2", width=0.5);
948// isects = circle_line_intersection(r=r, cp=cp, line=line, bounded=[true,false]);
949// color("#f44") move_copies(isects) circle(d=1);
950// Example(2D): If they don't intersect at all, then an empty list is returned.
951// line = [[-12,12], [12,8]];
952// cp = [-5,-2]; r = 10;
953// translate(cp) circle(r=r);
954// color("black") stroke(line, endcaps="arrow2", width=0.5);
955// isects = circle_line_intersection(r=r, cp=cp, line=line);
956// color("#f44") move_copies(isects) circle(d=1);
957function circle_line_intersection(r, cp, line, bounded=false, d, eps=EPSILON) =
958 assert(_valid_line(line,2), "Invalid 2d line.")
959 assert(is_vector(cp,2), "Circle center must be a 2-vector")
960 _circle_or_sphere_line_intersection(r, cp, line, bounded, d, eps);
961
962
963
964function _circle_or_sphere_line_intersection(r, cp, line, bounded=false, d, eps=EPSILON) =
965 let(r=get_radius(r=r,d=d,dflt=undef))
966 assert(is_num(r) && r>0, "Radius must be positive")
967 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition")
968 let(
969 bounded = force_list(bounded,2),
970 closest = line_closest_point(line,cp),
971 d = norm(closest-cp)
972 )
973 d > r ? [] :
974 let(
975 isect = approx(d,r,eps) ? [closest] :
976 let( offset = sqrt(r*r-d*d),
977 uvec=unit(line[1]-line[0])
978 ) [closest-offset*uvec, closest+offset*uvec]
979 )
980 [for(p=isect)
981 if ((!bounded[0] || (p-line[0])*(line[1]-line[0])>=0)
982 && (!bounded[1] || (p-line[1])*(line[0]-line[1])>=0)) p];
983
984
985// Function: circle_circle_intersection()
986// Synopsis: Find the intersection points of two 2d circles.
987// Topics: Geometry, Circles
988// See Also: circle_line_intersection(), circle_circle_intersection(), circle_2tangents(), circle_3points(), circle_point_tangents(), circle_circle_tangents()
989// Usage:
990// pts = circle_circle_intersection(r1|d1=, cp1, r2|d2=, cp2, [eps]);
991// Description:
992// Compute the intersection points of two circles. Returns a list of the intersection points, which
993// will contain two points in the general case, one point for tangent circles, or will be empty
994// if the circles do not intersect.
995// Arguments:
996// r1 = Radius of the first circle.
997// cp1 = Centerpoint of the first circle.
998// r2 = Radius of the second circle.
999// cp2 = Centerpoint of the second circle.
1000// eps = Tolerance for detecting tangent circles. Default: EPSILON
1001// ---
1002// d1 = Diameter of the first circle.
1003// d2 = Diameter of the second circle.
1004// Example(2D,NoAxes): Circles intersect in two points.
1005// $fn=32;
1006// cp1 = [4,4]; r1 = 3;
1007// cp2 = [7,7]; r2 = 2;
1008// pts = circle_circle_intersection(r1, cp1, r2, cp2);
1009// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1010// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1011// color("red") move_copies(pts) circle(r=.3);
1012// Example(2D,NoAxes): Circles are tangent, so one intersection point:
1013// $fn=32;
1014// cp1 = [4,4]; r1 = 4;
1015// cp2 = [4,10]; r2 = 2;
1016// pts = circle_circle_intersection(r1, cp1, r2, cp2);
1017// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1018// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1019// color("red") move_copies(pts) circle(r=.3);
1020// Example(2D,NoAxes): Another tangent example:
1021// $fn=32;
1022// cp1 = [4,4]; r1 = 4;
1023// cp2 = [5,5]; r2 = 4-sqrt(2);
1024// pts = circle_circle_intersection(r1, cp1, r2, cp2);
1025// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1026// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1027// color("red") move_copies(pts) circle(r=.3);
1028// Example(2D,NoAxes): Circles do not intersect. Returns empty list.
1029// $fn=32;
1030// cp1 = [3,4]; r1 = 2;
1031// cp2 = [7,10]; r2 = 3;
1032// pts = circle_circle_intersection(r1, cp1, r2, cp2);
1033// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1034// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1035// color("red") move_copies(pts) circle(r=.3);
1036function circle_circle_intersection(r1, cp1, r2, cp2, eps=EPSILON, d1, d2) =
1037 assert( is_path([cp1,cp2],dim=2), "Invalid center point(s)." )
1038 let(
1039 r1 = get_radius(r1=r1,d1=d1),
1040 r2 = get_radius(r1=r2,d1=d2),
1041 d = norm(cp2-cp1),
1042 a = (cp2-cp1)/d,
1043 b = [-a.y,a.x],
1044 L = (r1^2-r2^2+d^2)/2/d,
1045 hsqr = r1^2-L^2
1046 )
1047 approx(hsqr,0,eps) ? [L*a+cp1]
1048 : hsqr<0 ? []
1049 : let(h=sqrt(hsqr))
1050 [L*a+h*b+cp1, L*a-h*b+cp1];
1051
1052
1053// Function: circle_2tangents()
1054// Synopsis: Given two 2d or 3d rays, find a circle tangent to both.
1055// Topics: Geometry, Circles, Tangents
1056// See Also: circle_line_intersection(), circle_circle_intersection(), circle_2tangents(), circle_3points(), circle_point_tangents(), circle_circle_tangents()
1057// Usage:
1058// circ = circle_2tangents(r|d=, pt1, pt2, pt3, [tangents=]);
1059// circ = circle_2tangents(r|d=, [PT1, PT2, PT3], [tangents=]);
1060// Description:
1061// Given a pair of 2d or 3d rays with a common origin, and a known circle radius/diameter, finds
1062// the centerpoint for the circle of that size that touches both rays tangentally.
1063// Both rays start at `pt2`, one passing through `pt1`, and the other through `pt3`.
1064// .
1065// When called with collinear rays, returns `undef`.
1066// Otherwise, when called with `tangents=false`, returns `[CP,NORMAL]`.
1067// Otherwise, when called with `tangents=true`, returns `[CP,NORMAL,TANPT1,TANPT2]`.
1068// - CP is the centerpoint of the circle.
1069// - NORMAL is the normal vector of the plane that the circle is on (UP or DOWN if the points are 2D).
1070// - TANPT1 is the point where the circle is tangent to the ray `[pt2,pt1]`.
1071// - TANPT2 is the point where the circle is tangent to the ray `[pt2,pt3]`.
1072// Figure(3D,Med,NoAxes,VPD=130,VPT=[29,19,3],VPR=[55,0,25]):
1073// pts = [[45,10,-5], [10,5,10], [15,40,5]];
1074// rad = 15;
1075// circ = circle_2tangents(r=rad, pt1=pts[0], pt2=pts[1], pt3=pts[2], tangents=true);
1076// cp = circ[0]; n = circ[1]; tp1 = circ[2]; tp2 = circ[3];
1077// color("yellow") stroke(pts, endcaps="arrow2");
1078// color("purple") move_copies([cp,tp1,tp2]) sphere(d=2, $fn=12);
1079// color("lightgray") stroke([cp,tp2], width=0.5);
1080// stroke([cp,cp+n*20], endcap2="arrow2");
1081// labels = [
1082// ["pt1", "blue", 2.5, [ 4, 0, 1], pts[0]],
1083// ["pt2", "blue", 2.5, [-4, 0,-3], pts[1]],
1084// ["pt3", "blue", 2.5, [ 4, 0, 1], pts[2]],
1085// ["r", "blue", 2.5, [ 0,-2, 2], (cp+tp2)/2],
1086// ["CP", "brown", 2.5, [ 6,-4, 3], cp],
1087// ["Normal", "brown", 2.0, [ 5, 2, 1], cp+20*n],
1088// ["TanPt1", "brown", 2.0, [-5,-4, 0], tp1],
1089// ["TanPt2", "brown", 2.0, [-5, 0, 2], tp2],
1090// ];
1091// for(l=labels)
1092// color(l[1]) move(l[4]+l[3]) rot([55,0,25])
1093// linear_extrude(height=0.1)
1094// text(text=l[0], size=l[2], halign="center", valign="center");
1095// color("green",0.5) move(cp) cyl(h=0.1, r=rad, orient=n, $fn=36);
1096// Arguments:
1097// r = The radius of the circle to find.
1098// pt1 = A point that the first ray passes though.
1099// pt2 = The starting point of both rays.
1100// pt3 = A point that the second ray passes though.
1101// ---
1102// d = The diameter of the circle to find.
1103// tangents = If true, extended information about the tangent points is calculated and returned. Default: false
1104// Example(2D):
1105// pts = [[40,40], [10,10], [55,5]]; rad = 10;
1106// circ = circle_2tangents(r=rad, pt1=pts[0], pt2=pts[1], pt3=pts[2]);
1107// stroke(pts, endcaps="arrow2");
1108// color("red") move(circ[0]) circle(r=rad);
1109// Example(2D):
1110// pts = [[20,40], [10,10], [55,20]]; rad = 10;
1111// circ = circle_2tangents(r=rad, pt1=pts[0], pt2=pts[1], pt3=pts[2], tangents=true);
1112// stroke(pts, endcaps="arrow2");
1113// color("red") move(circ[0]) circle(r=rad);
1114// color("blue") move_copies(select(circ,2,3)) circle(d=2);
1115// Example(3D): Fit into 3D path corner.
1116// pts = [[45,5,10], [10,10,15], [30,40,30]]; rad = 10;
1117// circ = circle_2tangents(rad, [pts[0], pts[1], pts[2]]);
1118// stroke(pts, endcaps="arrow2");
1119// color("red") move(circ[0]) cyl(h=10, r=rad, orient=circ[1]);
1120// Example(3D):
1121// path = yrot(20, p=path3d(star(d=100, n=5, step=2)));
1122// stroke(path, closed=true);
1123// for (i = [0:1:5]) {
1124// crn = select(path, i*2-1, i*2+1);
1125// ci = circle_2tangents(5, crn[0], crn[1], crn[2]);
1126// move(ci[0]) cyl(h=10,r=5,orient=ci[1]);
1127// }
1128function circle_2tangents(r, pt1, pt2, pt3, tangents=false, d) =
1129 let(r = get_radius(r=r, d=d, dflt=undef))
1130 assert(r!=undef, "Must specify either r or d.")
1131 assert( ( is_path(pt1) && len(pt1)==3 && is_undef(pt2) && is_undef(pt3))
1132 || (is_matrix([pt1,pt2,pt3]) && (len(pt1)==2 || len(pt1)==3) ),
1133 str("Invalid input points. pt1=",pt1,", pt2=",pt2,", pt3=",pt3))
1134 is_undef(pt2)
1135 ? circle_2tangents(r, pt1[0], pt1[1], pt1[2], tangents=tangents)
1136 : is_collinear(pt1, pt2, pt3)? undef :
1137 let(
1138 v1 = unit(pt1 - pt2),
1139 v2 = unit(pt3 - pt2),
1140 vmid = unit(mean([v1, v2])),
1141 n = vector_axis(v1, v2),
1142 a = vector_angle(v1, v2),
1143 hyp = r / sin(a/2),
1144 cp = pt2 + hyp * vmid
1145 )
1146 !tangents ? [cp, n] :
1147 let(
1148 x = hyp * cos(a/2),
1149 tp1 = pt2 + x * v1,
1150 tp2 = pt2 + x * v2
1151 )
1152 [cp, n, tp1, tp2];
1153
1154
1155// Function: circle_3points()
1156// Synopsis: Find a circle passing through three 2d or 3d points.
1157// Topics: Geometry, Circles
1158// See Also: circle_line_intersection(), circle_circle_intersection(), circle_2tangents(), circle_3points(), circle_point_tangents(), circle_circle_tangents()
1159// Usage:
1160// circ = circle_3points(pt1, pt2, pt3);
1161// circ = circle_3points([PT1, PT2, PT3]);
1162// Description:
1163// Returns the [CENTERPOINT, RADIUS, NORMAL] of the circle that passes through three non-collinear
1164// points where NORMAL is the normal vector of the plane that the circle is on (UP or DOWN if the points are 2D).
1165// The centerpoint will be a 2D or 3D vector, depending on the points input. If all three
1166// points are 2D, then the resulting centerpoint will be 2D, and the normal will be UP ([0,0,1]).
1167// If any of the points are 3D, then the resulting centerpoint will be 3D. If the three points are
1168// collinear, then `[undef,undef,undef]` will be returned. The normal will be a normalized 3D
1169// vector with a non-negative Z axis. Instead of 3 arguments, it is acceptable to input the 3 points
1170// as a list given in `pt1`, leaving `pt2`and `pt3` as undef.
1171// Arguments:
1172// pt1 = The first point.
1173// pt2 = The second point.
1174// pt3 = The third point.
1175// Example(2D):
1176// pts = [[60,40], [10,10], [65,5]];
1177// circ = circle_3points(pts[0], pts[1], pts[2]);
1178// translate(circ[0]) color("green") stroke(circle(r=circ[1]),closed=true,$fn=72);
1179// translate(circ[0]) color("red") circle(d=3, $fn=12);
1180// move_copies(pts) color("blue") circle(d=3, $fn=12);
1181function circle_3points(pt1, pt2, pt3) =
1182 (is_undef(pt2) && is_undef(pt3) && is_list(pt1))
1183 ? circle_3points(pt1[0], pt1[1], pt1[2])
1184 : assert( is_vector(pt1) && is_vector(pt2) && is_vector(pt3)
1185 && max(len(pt1),len(pt2),len(pt3))<=3 && min(len(pt1),len(pt2),len(pt3))>=2,
1186 "Invalid point(s)." )
1187 is_collinear(pt1,pt2,pt3)? [undef,undef,undef] :
1188 let(
1189 v = [ point3d(pt1), point3d(pt2), point3d(pt3) ], // triangle vertices
1190 ed = [for(i=[0:2]) v[(i+1)%3]-v[i] ], // triangle edge vectors
1191 pm = [for(i=[0:2]) v[(i+1)%3]+v[i] ]/2, // edge mean points
1192 es = sortidx( [for(di=ed) norm(di) ] ),
1193 e1 = ed[es[1]], // take the 2 longest edges
1194 e2 = ed[es[2]],
1195 n0 = vector_axis(e1,e2), // normal standardization
1196 n = n0.z<0? -n0 : n0,
1197 sc = plane_intersection(
1198 [ each e1, e1*pm[es[1]] ], // planes orthogonal to 2 edges
1199 [ each e2, e2*pm[es[2]] ],
1200 [ each n, n*v[0] ]
1201 ), // triangle plane
1202 cp = len(pt1)+len(pt2)+len(pt3)>6 ? sc : [sc.x, sc.y],
1203 r = norm(sc-v[0])
1204 ) [ cp, r, n ];
1205
1206
1207
1208// Function: circle_point_tangents()
1209// Synopsis: Given a circle and point, find tangents to circle passing through the point.
1210// Topics: Geometry, Circles, Tangents
1211// See Also: circle_line_intersection(), circle_circle_intersection(), circle_2tangents(), circle_3points(), circle_point_tangents(), circle_circle_tangents()
1212// Usage:
1213// tangents = circle_point_tangents(r|d=, cp, pt);
1214// Description:
1215// Given a 2d circle and a 2d point outside that circle, finds the 2d tangent point(s) on the circle for a
1216// line passing through the point. Returns a list of zero or more 2D tangent points.
1217// Arguments:
1218// r = Radius of the circle.
1219// cp = The coordinates of the 2d circle centerpoint.
1220// pt = The coordinates of the 2d external point.
1221// ---
1222// d = Diameter of the circle.
1223// Example(2D):
1224// cp = [-10,-10]; r = 30; pt = [30,10];
1225// tanpts = circle_point_tangents(r=r, cp=cp, pt=pt);
1226// color("yellow") translate(cp) circle(r=r);
1227// color("cyan") for(tp=tanpts) {stroke([tp,pt]); stroke([tp,cp]);}
1228// color("red") move_copies(tanpts) circle(d=3,$fn=12);
1229// color("blue") move_copies([cp,pt]) circle(d=3,$fn=12);
1230function circle_point_tangents(r, cp, pt, d) =
1231 assert(is_finite(r) || is_finite(d), "Invalid radius or diameter." )
1232 assert(is_path([cp, pt],dim=2), "Invalid center point or external point.")
1233 let(
1234 r = get_radius(r=r, d=d, dflt=1),
1235 delta = pt - cp,
1236 dist = norm(delta),
1237 baseang = atan2(delta.y,delta.x)
1238 ) dist < r? [] :
1239 approx(dist,r)? [pt] :
1240 let(
1241 relang = acos(r/dist),
1242 angs = [baseang + relang, baseang - relang]
1243 ) [for (ang=angs) cp + r*[cos(ang),sin(ang)]];
1244
1245
1246// Function: circle_circle_tangents()
1247// Synopsis: Find tangents to a pair of circles in 2d.
1248// Topics: Geometry, Circles, Tangents
1249// See Also: circle_line_intersection(), circle_circle_intersection(), circle_2tangents(), circle_3points(), circle_point_tangents(), circle_circle_tangents()
1250// Usage:
1251// segs = circle_circle_tangents(r1|d1=, cp1, r2|d2=, cp2);
1252// Description:
1253// Computes 2d lines tangents to a pair of circles in 2d. Returns a list of line endpoints [p1,p2] where
1254// p1 is the tangent point on circle 1 and p2 is the tangent point on circle 2.
1255// If four tangents exist then the first one is the left hand exterior tangent as regarded looking from
1256// circle 1 toward circle 2. The second value is the right hand exterior tangent. The third entry
1257// gives the interior tangent that starts on the left of circle 1 and crosses to the right side of
1258// circle 2. And the fourth entry is the last interior tangent that starts on the right side of
1259// circle 1. If the circles intersect then the interior tangents don't exist and the function
1260// returns only two entries. If one circle is inside the other one then no tangents exist
1261// so the function returns the empty set. When the circles are tangent a degenerate tangent line
1262// passes through the point of tangency of the two circles: this degenerate line is NOT returned.
1263// Arguments:
1264// r1 = Radius of the first circle.
1265// cp1 = Centerpoint of the first circle.
1266// r2 = Radius of the second circle.
1267// cp2 = Centerpoint of the second circle.
1268// ---
1269// d1 = Diameter of the first circle.
1270// d2 = Diameter of the second circle.
1271// Example(2D,NoAxes): Four tangents, first in green, second in black, third in blue, last in red.
1272// $fn=32;
1273// cp1 = [3,4]; r1 = 2;
1274// cp2 = [7,10]; r2 = 3;
1275// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1276// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1277// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1278// colors = ["green","black","blue","red"];
1279// for(i=[0:len(pts)-1]) color(colors[i]) stroke(pts[i],width=0.2);
1280// Example(2D,NoAxes): Circles overlap so only exterior tangents exist.
1281// $fn=32;
1282// cp1 = [4,4]; r1 = 3;
1283// cp2 = [7,7]; r2 = 2;
1284// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1285// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1286// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1287// colors = ["green","black","blue","red"];
1288// for(i=[0:len(pts)-1]) color(colors[i]) stroke(pts[i],width=0.2);
1289// Example(2D,NoAxes): Circles are tangent. Only exterior tangents are returned. The degenerate internal tangent is not returned.
1290// $fn=32;
1291// cp1 = [4,4]; r1 = 4;
1292// cp2 = [4,10]; r2 = 2;
1293// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1294// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1295// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1296// colors = ["green","black","blue","red"];
1297// for(i=[0:1:len(pts)-1]) color(colors[i]) stroke(pts[i],width=0.2);
1298// Example(2D,NoAxes): One circle is inside the other: no tangents exist. If the interior circle is tangent the single degenerate tangent will not be returned.
1299// $fn=32;
1300// cp1 = [4,4]; r1 = 4;
1301// cp2 = [5,5]; r2 = 2;
1302// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1303// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1304// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1305// echo(pts); // Returns []
1306function circle_circle_tangents(r1, cp1, r2, cp2, d1, d2) =
1307 assert( is_path([cp1,cp2],dim=2), "Invalid center point(s)." )
1308 let(
1309 r1 = get_radius(r1=r1,d1=d1),
1310 r2 = get_radius(r1=r2,d1=d2),
1311 Rvals = [r2-r1, r2-r1, -r2-r1, -r2-r1]/norm(cp1-cp2),
1312 kvals = [-1,1,-1,1],
1313 ext = [1,1,-1,-1],
1314 N = 1-sqr(Rvals[2])>=0 ? 4 :
1315 1-sqr(Rvals[0])>=0 ? 2 : 0,
1316 coef= [
1317 for(i=[0:1:N-1]) [
1318 [Rvals[i], -kvals[i]*sqrt(1-sqr(Rvals[i]))],
1319 [kvals[i]*sqrt(1-sqr(Rvals[i])), Rvals[i]]
1320 ] * unit(cp2-cp1)
1321 ]
1322 ) [
1323 for(i=[0:1:N-1]) let(
1324 pt = [
1325 cp1-r1*coef[i],
1326 cp2-ext[i]*r2*coef[i]
1327 ]
1328 ) if (pt[0]!=pt[1]) pt
1329 ];
1330
1331
1332
1333/// Internal Function: _noncollinear_triple()
1334/// Usage:
1335/// bool = _noncollinear_triple(points);
1336/// Topics: Geometry, Noncollinearity
1337/// Description:
1338/// Finds the indices of three non-collinear points from the pointlist `points`.
1339/// It selects two well separated points to define a line and chooses the third point
1340/// to be the point farthest off the line. The points do not necessarily having the
1341/// same winding direction as the polygon so they cannot be used to determine the
1342/// winding direction or the direction of the normal.
1343/// If all points are collinear returns [] when `error=true` or an error otherwise .
1344/// Arguments:
1345/// points = List of input points.
1346/// error = Defines the behaviour for collinear input points. When `true`, produces an error, otherwise returns []. Default: `true`.
1347/// eps = Tolerance for collinearity test. Default: EPSILON.
1348function _noncollinear_triple(points,error=true,eps=EPSILON) =
1349 assert( is_path(points), "Invalid input points." )
1350 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
1351 len(points)<3 ? [] :
1352 let(
1353 pa = points[0],
1354 b = furthest_point(pa, points),
1355 pb = points[b],
1356 nrm = norm(pa-pb)
1357 )
1358 nrm <= eps ?
1359 assert(!error, "Cannot find three noncollinear points in pointlist.") [] :
1360 let(
1361 n = (pb-pa)/nrm,
1362 distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)]
1363 )
1364 max(distlist) < eps*nrm ?
1365 assert(!error, "Cannot find three noncollinear points in pointlist.") [] :
1366 [0, b, max_index(distlist)];
1367
1368
1369
1370// Section: Sphere Calculations
1371
1372
1373// Function: sphere_line_intersection()
1374// Synopsis: Find intersection between a sphere and line, ray or segment.
1375// Topics: Geometry, Spheres, Lines, Intersection
1376// See Also: circle_line_intersection(), circle_circle_intersection(), circle_2tangents(), circle_3points(), circle_point_tangents(), circle_circle_tangents()
1377// Usage:
1378// isect = sphere_line_intersection(r|d=, cp, line, [bounded], [eps=]);
1379// Description:
1380// Find intersection points between a sphere and a line, ray or segment specified by two points.
1381// By default the line is unbounded.
1382// Arguments:
1383// r = Radius of sphere
1384// cp = Centerpoint of sphere
1385// line = Two points defining the line
1386// bounded = false for unbounded line, true for a segment, or a vector [false,true] or [true,false] to specify a ray with the first or second end unbounded. Default: false
1387// ---
1388// d = diameter of sphere
1389// eps = epsilon used for identifying the case with one solution. Default: 1e-9
1390// Example(3D):
1391// cp = [10,20,5]; r = 40;
1392// line = [[-50,-10,25], [70,0,40]];
1393// isects = sphere_line_intersection(r=r, cp=cp, line=line);
1394// color("cyan") stroke(line);
1395// move(cp) sphere(r=r, $fn=72);
1396// color("red") move_copies(isects) sphere(d=3, $fn=12);
1397function sphere_line_intersection(r, cp, line, bounded=false, d, eps=EPSILON) =
1398 assert(_valid_line(line,3), "Invalid 3d line.")
1399 assert(is_vector(cp,3), "Sphere center must be a 3-vector")
1400 _circle_or_sphere_line_intersection(r, cp, line, bounded, d, eps);
1401
1402
1403
1404
1405// Section: Polygons
1406
1407// Function: polygon_area()
1408// Synopsis: Calculate area of a 2d or 3d polygon.
1409// Topics: Geometry, Polygons, Area
1410// See Also: polygon_area(), centroid(), polygon_normal(), point_in_polygon(), polygon_line_intersection()
1411// Usage:
1412// area = polygon_area(poly, [signed]);
1413// Description:
1414// Given a 2D or 3D simple planar polygon, returns the area of that polygon.
1415// If the polygon is non-planar the result is `undef.` If the polygon is self-intersecting
1416// then the return will be a meaningless number.
1417// When `signed` is true and the polygon is 2d, a signed area is returned: a positive area indicates a counter-clockwise polygon.
1418// The area of 3d polygons is always nonnegative.
1419// Arguments:
1420// poly = Polygon to compute the area of.
1421// signed = If true, a signed area is returned. Default: false.
1422function polygon_area(poly, signed=false) =
1423 assert(is_path(poly), "Invalid polygon." )
1424 len(poly)<3 ? 0 :
1425 len(poly)==3 ?
1426 let( total= len(poly[0])==2 ? 0.5*cross(poly[2]-poly[0],poly[2]-poly[1]) : 0.5*norm(cross(poly[2]-poly[0],poly[2]-poly[1])))
1427 signed ? total : abs(total) :
1428 len(poly[0])==2
1429 ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 )
1430 signed ? total : abs(total)
1431 : let( plane = plane_from_polygon(poly) )
1432 is_undef(plane) ? undef :
1433 let(
1434 n = plane_normal(plane),
1435 total =
1436 -sum([ for(i=[1:1:len(poly)-2])
1437 cross(poly[i]-poly[0], poly[i+1]-poly[0])
1438 ]) * n/2
1439 )
1440 signed ? total : abs(total);
1441
1442
1443// Function: centroid()
1444// Synopsis: Compute centroid of a 2d or 3d polygon or a VNF.
1445// Topics: Geometry, Polygons, Centroid
1446// See Also: polygon_area(), centroid(), polygon_normal(), point_in_polygon(), polygon_line_intersection()
1447// Usage:
1448// c = centroid(object, [eps]);
1449// Description:
1450// Given a simple 2D polygon, returns the 2D coordinates of the polygon's centroid.
1451// Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
1452// If you provide a non-planar or collinear polygon you will get an error. For self-intersecting
1453// polygons you may get an error or you may get meaningless results.
1454// .
1455// Given a [region](regions.scad), returns the 2D coordinates of the region's centroid.
1456// .
1457// Given a manifold [VNF](vnf.scad) then returns the 3D centroid of the polyhedron. The VNF must
1458// describe a valid polyhedron with consistent face direction and no holes in the mesh; otherwise
1459// the results are undefined.
1460// Arguments:
1461// object = object to compute the centroid of
1462// eps = epsilon value for identifying degenerate cases
1463// Example(2D):
1464// path = [
1465// [-10,10], [-5,15], [15,15], [20,0],
1466// [15,-5], [25,-20], [25,-27], [15,-20],
1467// [0,-30], [-15,-25], [-5,-5]
1468// ];
1469// linear_extrude(height=0.01) polygon(path);
1470// cp = centroid(path);
1471// color("red") move(cp) sphere(d=2);
1472function centroid(object,eps=EPSILON) =
1473 assert(is_finite(eps) && (eps>=0), "The tolerance should a non-negative value." )
1474 is_vnf(object) ? _vnf_centroid(object,eps)
1475 : is_path(object,[2,3]) ? _polygon_centroid(object,eps)
1476 : is_region(object) ? (len(object)==1 ? _polygon_centroid(object[0],eps) : _region_centroid(object,eps))
1477 : assert(false, "Input must be a VNF, a region, or a 2D or 3D polygon");
1478
1479
1480/// Internal Function: _region_centroid()
1481/// Compute centroid of region
1482function _region_centroid(region,eps=EPSILON) =
1483 let(
1484 region=force_region(region),
1485 parts = region_parts(region),
1486 // Rely on region_parts returning all outside polygons clockwise
1487 // and inside (hole) polygons counterclockwise, so areas have reversed sign
1488 cent_area = [for(R=parts, p=R)
1489 let(A=polygon_area(p,signed=true))
1490 [A*_polygon_centroid(p),A]],
1491 total = sum(cent_area)
1492 )
1493 total[0]/total[1];
1494
1495
1496/// Internal Function: _polygon_centroid()
1497/// Usage:
1498/// cpt = _polygon_centroid(poly);
1499/// Topics: Geometry, Polygons, Centroid
1500/// Description:
1501/// Given a simple 2D polygon, returns the 2D coordinates of the polygon's centroid.
1502/// Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
1503/// Collinear points produce an error. The results are meaningless for self-intersecting
1504/// polygons or an error is produced.
1505/// Arguments:
1506/// poly = Points of the polygon from which the centroid is calculated.
1507/// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
1508function _polygon_centroid(poly, eps=EPSILON) =
1509 assert( is_path(poly,dim=[2,3]), "The input must be a 2D or 3D polygon." )
1510 let(
1511 n = len(poly[0])==2 ? 1 :
1512 let( plane = plane_from_points(poly, fast=false))
1513 assert(!is_undef(plane), "The polygon must be planar." )
1514 plane_normal(plane),
1515 v0 = poly[0] ,
1516 val = sum([
1517 for(i=[1:len(poly)-2])
1518 let(
1519 v1 = poly[i],
1520 v2 = poly[i+1],
1521 area = cross(v2-v0,v1-v0)*n
1522 ) [ area, (v0+v1+v2)*area ]
1523 ])
1524 )
1525 assert(!approx(val[0],0, eps), "The polygon is self-intersecting or its points are collinear.")
1526 val[1]/val[0]/3;
1527
1528
1529
1530// Function: polygon_normal()
1531// Synopsis: Return normal to a polygon.
1532// Topics: Geometry, Polygons
1533// See Also: polygon_area(), centroid(), polygon_normal(), point_in_polygon(), polygon_line_intersection()
1534// Usage:
1535// vec = polygon_normal(poly);
1536// Description:
1537// Given a 3D simple planar polygon, returns a unit normal vector for the polygon. The vector
1538// is oriented so that if the normal points towards the viewer, the polygon winds in the clockwise
1539// direction. If the polygon has zero area, returns `undef`. If the polygon is self-intersecting
1540// the the result is undefined. It doesn't check for coplanarity.
1541// Arguments:
1542// poly = The list of 3D path points for the perimeter of the polygon.
1543// Example(3D):
1544// path = rot([0,30,15], p=path3d(star(n=5, d=100, step=2)));
1545// stroke(path, closed=true);
1546// n = polygon_normal(path);
1547// rot(from=UP, to=n)
1548// color("red")
1549// stroke([[0,0,0], [0,0,20]], endcap2="arrow2");
1550function polygon_normal(poly) =
1551 assert(is_path(poly,dim=3), "Invalid 3D polygon." )
1552 let(
1553 area_vec = sum([for(i=[1:len(poly)-2])
1554 cross(poly[i]-poly[0],
1555 poly[i+1]-poly[i])])
1556 )
1557 unit(-area_vec, error=undef);
1558
1559
1560// Function: point_in_polygon()
1561// Synopsis: Checks if a 2d point is inside or on the boundary of a 2d polygon.
1562// Topics: Geometry, Polygons
1563// See Also: polygon_area(), centroid(), polygon_normal(), point_in_polygon(), polygon_line_intersection()
1564// Usage:
1565// bool = point_in_polygon(point, poly, [nonzero], [eps])
1566// Description:
1567// This function tests whether the given 2D point is inside, outside or on the boundary of
1568// the specified 2D polygon.
1569// The polygon is given as a list of 2D points, not including the repeated end point.
1570// Returns -1 if the point is outside the polygon.
1571// Returns 0 if the point is on the boundary.
1572// Returns 1 if the point lies in the interior.
1573// The polygon does not need to be simple: it may have self-intersections.
1574// But the polygon cannot have holes (it must be simply connected).
1575// Rounding errors may give mixed results for points on or near the boundary.
1576// .
1577// When polygons intersect themselves different definitions exist for determining which points
1578// are inside the polygon. The figure below shows the difference.
1579// OpenSCAD uses the Even-Odd rule when creating polygons, where membership in overlapping regions
1580// depends on how many times they overlap. The Nonzero rule considers point inside the polygon if
1581// the polygon overlaps them any number of times. For more information see
1582// https://en.wikipedia.org/wiki/Nonzero-rule and https://en.wikipedia.org/wiki/Even–odd_rule.
1583// Figure(2D,Med,NoAxes):
1584// a=20;
1585// b=30;
1586// ofs = 17;
1587// curve = [for(theta=[0:10:140]) [a * theta/360*2*PI - b*sin(theta), a-b*cos(theta)-20]];
1588// path = deduplicate(concat( reverse(offset(curve,r=ofs)),
1589// xflip(offset(curve,r=ofs)),
1590// xflip(reverse(curve)),
1591// curve
1592// ));
1593// left(40){
1594// polygon(path);
1595// color("red")stroke(path, width=1, closed=true);
1596// color("red")back(28/(2/3))text("Even-Odd", size=5/(2/3), halign="center");
1597// }
1598// right(40){
1599// dp = polygon_parts(path,nonzero=true);
1600// region(dp);
1601// color("red"){stroke(path,width=1,closed=true);
1602// back(28/(2/3))text("Nonzero", size=5/(2/3), halign="center");
1603// }
1604// }
1605// Arguments:
1606// point = The 2D point to check
1607// poly = The list of 2D points forming the perimeter of the polygon.
1608// nonzero = The rule to use: true for "Nonzero" rule and false for "Even-Odd". Default: false (Even-Odd)
1609// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
1610// Example(2D): With nonzero set to false (the default), we get this result. Green dots are inside the polygon and red are outside:
1611// a=20*2/3;
1612// b=30*2/3;
1613// ofs = 17*2/3;
1614// curve = [for(theta=[0:10:140]) [a * theta/360*2*PI - b*sin(theta), a-b*cos(theta)]];
1615// path = deduplicate(concat( reverse(offset(curve,r=ofs)),
1616// xflip(offset(curve,r=ofs)),
1617// xflip(reverse(curve)),
1618// curve
1619// ));
1620// stroke(path,closed=true);
1621// pts = [[0,0],[10,0],[0,20]];
1622// for(p=pts){
1623// color(point_in_polygon(p,path)==1 ? "green" : "red")
1624// move(p)circle(r=1.5, $fn=12);
1625// }
1626// Example(2D): With nonzero set to true, one dot changes color:
1627// a=20*2/3;
1628// b=30*2/3;
1629// ofs = 17*2/3;
1630// curve = [for(theta=[0:10:140]) [a * theta/360*2*PI - b*sin(theta), a-b*cos(theta)]];
1631// path = deduplicate(concat( reverse(offset(curve,r=ofs)),
1632// xflip(offset(curve,r=ofs)),
1633// xflip(reverse(curve)),
1634// curve
1635// ));
1636// stroke(path,closed=true);
1637// pts = [[0,0],[10,0],[0,20]];
1638// for(p=pts){
1639// color(point_in_polygon(p,path,nonzero=true)==1 ? "green" : "red")
1640// move(p)circle(r=1.5, $fn=12);
1641// }
1642
1643// Internal function for point_in_polygon
1644
1645function _point_above_below_segment(point, edge) =
1646 let( edge = edge - [point, point] )
1647 edge[0].y <= 0
1648 ? (edge[1].y > 0 && cross(edge[0], edge[1]-edge[0]) > 0) ? 1 : 0
1649 : (edge[1].y <= 0 && cross(edge[0], edge[1]-edge[0]) < 0) ? -1 : 0;
1650
1651
1652function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
1653 // Original algorithms from http://geomalgorithms.com/a03-_inclusion.html
1654 assert( is_vector(point,2) && is_path(poly,dim=2) && len(poly)>2,
1655 "The point and polygon should be in 2D. The polygon should have more that 2 points." )
1656 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
1657 // Check bounding box
1658 let(
1659 box = pointlist_bounds(poly)
1660 )
1661 point.x<box[0].x-eps || point.x>box[1].x+eps
1662 || point.y<box[0].y-eps || point.y>box[1].y+eps ? -1
1663 :
1664 // Does the point lie on any edges? If so return 0.
1665 let(
1666 segs = pair(poly,true),
1667 on_border = [for (seg=segs)
1668 if (norm(seg[0]-seg[1])>eps && _is_point_on_line(point, seg, SEGMENT, eps=eps)) 1]
1669 )
1670 on_border != [] ? 0 :
1671 nonzero // Compute winding number and return 1 for interior, -1 for exterior
1672 ? let(
1673 winding = [
1674 for(seg=segs)
1675 let(
1676 p0=seg[0]-point,
1677 p1=seg[1]-point
1678 )
1679 if (norm(p0-p1)>eps)
1680 p0.y <=0
1681 ? p1.y > 0 && cross(p0,p1-p0)>0 ? 1 : 0
1682 : p1.y <=0 && cross(p0,p1-p0)<0 ? -1: 0
1683 ]
1684 )
1685 sum(winding) != 0 ? 1 : -1
1686 : // or compute the crossings with the ray [point, point+[1,0]]
1687 let(
1688 cross = [
1689 for(seg=segs)
1690 let(
1691 p0 = seg[0]-point,
1692 p1 = seg[1]-point
1693 )
1694 if (
1695 ( (p1.y>eps && p0.y<=eps) || (p1.y<=eps && p0.y>eps) )
1696 && -eps < p0.x - p0.y *(p1.x - p0.x)/(p1.y - p0.y)
1697 )
1698 1
1699 ]
1700 )
1701 2*(len(cross)%2)-1;
1702
1703
1704
1705// Function: polygon_line_intersection()
1706// Synopsis: Find intersection between 2d or 3d polygon and a line, segment or ray.
1707// Topics: Geometry, Polygons, Lines, Intersection
1708// See Also: polygon_area(), centroid(), polygon_normal(), point_in_polygon(), polygon_line_intersection()
1709// Usage:
1710// pt = polygon_line_intersection(poly, line, [bounded], [nonzero], [eps]);
1711// Description:
1712// Takes a possibly bounded line, and a 2D or 3D planar polygon, and finds their intersection. Note the polygon is
1713// treated as its boundary and interior, so the intersection may include both points and line segments.
1714// If the line does not intersect the polygon then returns `undef`.
1715// In 3D if the line is not on the plane of the polygon but intersects it then you get a single intersection point.
1716// Otherwise the polygon and line are in the same plane, or when your input is 2D, you will get a list of segments and
1717// single point lists. Use `is_vector` to distinguish these two cases.
1718// .
1719// In the 2D case, a common result is a list containing a single segment, which lists the two intersection points
1720// with the boundary of the polygon.
1721// When single points are in the intersection (the line just touches a polygon corner) they appear on the segment
1722// list as lists of a single point
1723// (like single point segments) so a single point intersection in 2D has the form `[[[x,y,z]]]` as compared
1724// to a single point intersection in 3D which has the form `[x,y,z]`. You can identify whether an entry in the
1725// segment list is a true segment by checking its length, which will be 2 for a segment and 1 for a point.
1726// Arguments:
1727// poly = The 3D planar polygon to find the intersection with.
1728// line = A list of two distinct 3D points on the line.
1729// bounded = If false, the line is considered unbounded. If true, it is treated as a bounded line segment. If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray. Default: false (unbounded)
1730// nonzero = set to true to use the nonzero rule for determining it points are in a polygon. See point_in_polygon. Default: false.
1731// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
1732// Example(3D): The line intersects the 3d hexagon in a single point.
1733// hex = zrot(140,p=rot([-45,40,20],p=path3d(hexagon(r=15))));
1734// line = [[5,0,-13],[-3,-5,13]];
1735// isect = polygon_line_intersection(hex,line);
1736// stroke(hex,closed=true);
1737// stroke(line);
1738// color("red")move(isect)sphere(r=1,$fn=12);
1739// Example(2D): In 2D things are more complicated. The output is a list of intersection parts, in the simplest case a single segment.
1740// hex = hexagon(r=15);
1741// line = [[-20,10],[25,-7]];
1742// isect = polygon_line_intersection(hex,line);
1743// stroke(hex,closed=true);
1744// stroke(line,endcaps="arrow2");
1745// color("red")
1746// for(part=isect)
1747// if(len(part)==1)
1748// move(part[0]) sphere(r=1);
1749// else
1750// stroke(part);
1751// Example(2D): Here the line is treated as a ray.
1752// hex = hexagon(r=15);
1753// line = [[0,0],[25,-7]];
1754// isect = polygon_line_intersection(hex,line,RAY);
1755// stroke(hex,closed=true);
1756// stroke(line,endcap2="arrow2");
1757// color("red")
1758// for(part=isect)
1759// if(len(part)==1)
1760// move(part[0]) circle(r=1,$fn=12);
1761// else
1762// stroke(part);
1763// Example(2D): Here the intersection is a single point, which is returned as a single point "path" on the path list.
1764// hex = hexagon(r=15);
1765// line = [[15,-10],[15,13]];
1766// isect = polygon_line_intersection(hex,line,RAY);
1767// stroke(hex,closed=true);
1768// stroke(line,endcap2="arrow2");
1769// color("red")
1770// for(part=isect)
1771// if(len(part)==1)
1772// move(part[0]) circle(r=1,$fn=12);
1773// else
1774// stroke(part);
1775// Example(2D): Another way to get a single segment
1776// hex = hexagon(r=15);
1777// line = rot(30,p=[[15,-10],[15,25]],cp=[15,0]);
1778// isect = polygon_line_intersection(hex,line,RAY);
1779// stroke(hex,closed=true);
1780// stroke(line,endcap2="arrow2");
1781// color("red")
1782// for(part=isect)
1783// if(len(part)==1)
1784// move(part[0]) circle(r=1,$fn=12);
1785// else
1786// stroke(part);
1787// Example(2D): Single segment again
1788// star = star(r=15,n=8,step=2);
1789// line = [[20,-5],[-5,20]];
1790// isect = polygon_line_intersection(star,line,RAY);
1791// stroke(star,closed=true);
1792// stroke(line,endcap2="arrow2");
1793// color("red")
1794// for(part=isect)
1795// if(len(part)==1)
1796// move(part[0]) circle(r=1,$fn=12);
1797// else
1798// stroke(part);
1799// Example(2D): Solution is two points
1800// star = star(r=15,n=8,step=3);
1801// line = rot(22.5,p=[[15,-10],[15,20]],cp=[15,0]);
1802// isect = polygon_line_intersection(star,line,SEGMENT);
1803// stroke(star,closed=true);
1804// stroke(line);
1805// color("red")
1806// for(part=isect)
1807// if(len(part)==1)
1808// move(part[0]) circle(r=1,$fn=12);
1809// else
1810// stroke(part);
1811// Example(2D): Solution is list of three segments
1812// star = star(r=25,ir=9,n=8);
1813// line = [[-25,12],[25,12]];
1814// isect = polygon_line_intersection(star,line);
1815// stroke(star,closed=true);
1816// stroke(line,endcaps="arrow2");
1817// color("red")
1818// for(part=isect)
1819// if(len(part)==1)
1820// move(part[0]) circle(r=1,$fn=12);
1821// else
1822// stroke(part);
1823// Example(2D): Solution is a mixture of segments and points
1824// star = star(r=25,ir=9,n=7);
1825// line = [left(10,p=star[8]), right(50,p=star[8])];
1826// isect = polygon_line_intersection(star,line);
1827// stroke(star,closed=true);
1828// stroke(line,endcaps="arrow2");
1829// color("red")
1830// for(part=isect)
1831// if(len(part)==1)
1832// move(part[0]) circle(r=1,$fn=12);
1833// else
1834// stroke(part);
1835function polygon_line_intersection(poly, line, bounded=false, nonzero=false, eps=EPSILON) =
1836 assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
1837 assert(is_path(poly,dim=[2,3]), "Invalid polygon." )
1838 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.")
1839 assert(_valid_line(line,dim=len(poly[0]),eps=eps), "Line invalid or does not match polygon dimension." )
1840 let(
1841 bounded = force_list(bounded,2),
1842 poly = deduplicate(poly)
1843 )
1844 len(poly[0])==2 ? // planar case
1845 let(
1846 linevec = unit(line[1] - line[0]),
1847 bound = 100*max(v_abs(flatten(pointlist_bounds(poly)))),
1848 boundedline = [line[0] + (bounded[0]? 0 : -bound) * linevec,
1849 line[1] + (bounded[1]? 0 : bound) * linevec],
1850 parts = split_region_at_region_crossings(boundedline, [poly], closed1=false)[0][0],
1851 inside = [
1852 if(point_in_polygon(parts[0][0], poly, nonzero=nonzero, eps=eps) == 0)
1853 [parts[0][0]], // Add starting point if it is on the polygon
1854 for(part = parts)
1855 if (point_in_polygon(mean(part), poly, nonzero=nonzero, eps=eps) >=0 )
1856 part
1857 else if(len(part)==2 && point_in_polygon(part[1], poly, nonzero=nonzero, eps=eps) == 0)
1858 [part[1]] // Add segment end if it is on the polygon
1859 ]
1860 )
1861 (len(inside)==0 ? undef : _merge_segments(inside, [inside[0]], eps))
1862 : // 3d case
1863 let(indices = _noncollinear_triple(poly))
1864 indices==[] ? undef : // Polygon is collinear
1865 let(
1866 plane = plane3pt(poly[indices[0]], poly[indices[1]], poly[indices[2]]),
1867 plane_isect = plane_line_intersection(plane, line, bounded, eps)
1868 )
1869 is_undef(plane_isect) ? undef :
1870 is_vector(plane_isect,3) ?
1871 let(
1872 poly2d = project_plane(plane,poly),
1873 pt2d = project_plane(plane, plane_isect)
1874 )
1875 (point_in_polygon(pt2d, poly2d, nonzero=nonzero, eps=eps) < 0 ? undef : plane_isect)
1876 : // Case where line is on the polygon plane
1877 let(
1878 poly2d = project_plane(plane, poly),
1879 line2d = project_plane(plane, line),
1880 segments = polygon_line_intersection(poly2d, line2d, bounded=bounded, nonzero=nonzero, eps=eps)
1881 )
1882 segments==undef ? undef
1883 : [for(seg=segments) len(seg)==2 ? lift_plane(plane,seg) : [lift_plane(plane,seg[0])]];
1884
1885function _merge_segments(insegs,outsegs, eps, i=1) =
1886 i==len(insegs) ? outsegs :
1887 approx(last(last(outsegs)), insegs[i][0], eps)
1888 ? _merge_segments(insegs, [each list_head(outsegs),[last(outsegs)[0],last(insegs[i])]], eps, i+1)
1889 : _merge_segments(insegs, [each outsegs, insegs[i]], eps, i+1);
1890
1891
1892
1893// Function: polygon_triangulate()
1894// Synopsis: Divide a polygon into triangles.
1895// Topics: Geometry, Triangulation
1896// See Also: vnf_triangulate()
1897// Usage:
1898// triangles = polygon_triangulate(poly, [ind], [error], [eps])
1899// Description:
1900// Given a simple polygon in 2D or 3D, triangulates it and returns a list
1901// of triples indexing into the polygon vertices. When the optional argument `ind` is
1902// given, it is used as an index list into `poly` to define the polygon vertices. In that case,
1903// `poly` may have a length greater than `ind`. When `ind` is undefined, all points in `poly`
1904// are considered as vertices of the polygon.
1905// .
1906// For 2d polygons, the output triangles will have the same winding (CW or CCW) of
1907// the input polygon. For 3d polygons, the triangle windings will induce a normal
1908// vector with the same direction of the polygon normal.
1909// .
1910// The function produce correct triangulations for some non-twisted non-simple polygons.
1911// A polygon is non-twisted iff it is simple or it has a partition in
1912// simple polygons with the same winding such that the intersection of any two partitions is
1913// made of full edges and/or vertices of both partitions. These polygons may have "touching" vertices
1914// (two vertices having the same coordinates, but distinct adjacencies) and "contact" edges
1915// (edges whose vertex pairs have the same pairwise coordinates but are in reversed order) but has
1916// no self-crossing. See examples bellow. If all polygon edges are contact edges (polygons with
1917// zero area), it returns an empty list for 2d polygons and reports an error for 3d polygons.
1918// Triangulation errors are reported either by an assert error (when `error=true`) or by returning
1919// `undef` (when `error=false`). Invalid arguments always produce an assert error.
1920// .
1921// Twisted polygons have no consistent winding and when input to this function usually reports
1922// an error but when an error is not reported the outputs are not correct triangulations. The function
1923// can work for 3d non-planar polygons if they are close enough to planar but may otherwise
1924// report an error for this case.
1925// Arguments:
1926// poly = Array of the polygon vertices.
1927// ind = If given, a list of indices indexing the vertices of the polygon in `poly`. Default: use all the points of poly
1928// error = If false, returns `undef` when the polygon cannot be triangulated; otherwise, issues an assert error. Default: true.
1929// eps = A maximum tolerance in geometrical tests. Default: EPSILON
1930// Example(2D,NoAxes): a simple polygon; see from above
1931// poly = star(id=10, od=15,n=11);
1932// tris = polygon_triangulate(poly);
1933// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1934// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1935// color("magenta") up(2) stroke(poly,.25,closed=true);
1936// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1937// Example(2D,NoAxes): a polygon with a hole and one "contact" edge; see from above
1938// poly = [ [-10,0], [10,0], [0,10], [-10,0], [-4,4], [4,4], [0,2], [-4,4] ];
1939// tris = polygon_triangulate(poly);
1940// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1941// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1942// color("magenta") up(2) stroke(poly,.25,closed=true);
1943// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1944// Example(2D,NoAxes): a polygon with "touching" vertices and no holes; see from above
1945// poly = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ];
1946// tris = polygon_triangulate(poly);
1947// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1948// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1949// color("magenta") up(2) stroke(poly,.25,closed=true);
1950// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1951// Example(2D,NoAxes): a polygon with "contact" edges and no holes; see from above
1952// poly = [ [0,0], [10,0], [10,10], [0,10], [0,0], [3,3], [7,3],
1953// [7,7], [7,3], [3,3] ];
1954// tris = polygon_triangulate(poly);
1955// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1956// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1957// color("magenta") up(2) stroke(poly,.25,closed=true);
1958// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1959// Example(3D):
1960// include <BOSL2/polyhedra.scad>
1961// vnf = regular_polyhedron_info(name="dodecahedron",side=5,info="vnf");
1962// vnf_polyhedron(vnf);
1963// vnf_tri = [vnf[0], [for(face=vnf[1]) each polygon_triangulate(vnf[0], face) ] ];
1964// color("blue")
1965// vnf_wireframe(vnf_tri, width=.15);
1966function polygon_triangulate(poly, ind, error=true, eps=EPSILON) =
1967 assert(is_path(poly) && len(poly)>=3, "Polygon `poly` should be a list of at least three 2d or 3d points")
1968 assert(is_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
1969 "Improper or out of bounds list of indices")
1970 let( ind = is_undef(ind) ? count(len(poly)) : ind )
1971 len(ind) <=2 ? [] :
1972 len(ind) == 3
1973 ? _degenerate_tri([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] :
1974 // non zero area
1975 let( degen = norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) < 2*eps )
1976 assert( ! error || ! degen, "The polygon vertices are collinear.")
1977 degen ? undef : [ind]
1978 : len(poly[ind[0]]) == 3
1979 ? // find a representation of the polygon as a 2d polygon by projecting it on its own plane
1980 let(
1981 ind = deduplicate_indexed(poly, ind, eps)
1982 )
1983 len(ind)<3 ? [] :
1984 let(
1985 pts = select(poly,ind),
1986 nrm = -polygon_normal(pts)
1987 )
1988 assert( ! error || (nrm != undef),
1989 "The polygon has self-intersections or zero area or its vertices are collinear or non coplanar.")
1990 nrm == undef ? undef :
1991 let(
1992 imax = max_index([for(p=pts) norm(p-pts[0]) ]),
1993 v1 = unit( pts[imax] - pts[0] ),
1994 v2 = cross(v1,nrm),
1995 prpts = pts*transpose([v1,v2]) // the 2d projection of pts on the polygon plane
1996 )
1997 let( tris = _triangulate(prpts, count(len(ind)), error, eps) )
1998 tris == undef ? undef :
1999 [for(tri=tris) select(ind,tri) ]
2000 : is_polygon_clockwise(select(poly, ind))
2001 ? _triangulate( poly, ind, error, eps )
2002 : let( tris = _triangulate( poly, reverse(ind), error, eps ) )
2003 tris == undef ? undef :
2004 [for(tri=tris) reverse(tri) ];
2005
2006
2007// poly is supposed to be a 2d cw polygon
2008// implements a modified version of ear cut method for non-twisted polygons
2009// the polygons accepted by this function are those decomposable in simple
2010// CW polygons.
2011function _triangulate(poly, ind, error, eps=EPSILON, tris=[]) =
2012 len(ind)==3
2013 ? _degenerate_tri(select(poly,ind),eps)
2014 ? tris // if last 3 pts perform a degenerate triangle, ignore it
2015 : concat(tris,[ind]) // otherwise, include it
2016 : let( ear = _get_ear(poly,ind,eps) )
2017 assert( ! error || (ear != undef),
2018 "The polygon has twists or all its vertices are collinear or non coplanar.")
2019 ear == undef ? undef :
2020 is_list(ear) // is it a degenerate ear ?
2021 ? len(ind) <= 4 ? tris :
2022 _triangulate(poly, select(ind,ear[0]+3, ear[0]), error, eps, tris) // discard it
2023 : let(
2024 ear_tri = select(ind,ear,ear+2),
2025 indr = select(ind,ear+2, ear) // indices of the remaining path
2026 )
2027 _triangulate(poly, indr, error, eps, concat(tris,[ear_tri]));
2028
2029
2030// a returned ear will be:
2031// 1. a CW non-reflex triangle, made of subsequent poly vertices, without any other
2032// poly points inside except possibly at its own vertices
2033// 2. or a degenerate triangle where two vertices are coincident
2034// the returned ear is specified by the index of `ind` of its first vertex
2035function _get_ear(poly, ind, eps, _i=0) =
2036 let( lind = len(ind) )
2037 lind==3 ? 0 :
2038 let( // the _i-th ear candidate
2039 p0 = poly[ind[_i]],
2040 p1 = poly[ind[(_i+1)%lind]],
2041 p2 = poly[ind[(_i+2)%lind]]
2042 )
2043 // if vertex p1 is a convex candidate to be an ear,
2044 // check if the triangle [p0,p1,p2] contains any other point
2045 // except possibly p0 and p2
2046 // exclude the ear candidate central vertex p1 from the verts to check
2047 _tri_class([p0,p1,p2],eps) > 0
2048 && _none_inside(select(ind,_i+2, _i),poly,p0,p1,p2,eps) ? _i : // found an ear
2049 // otherwise check the next ear candidate
2050 _i<lind-1 ? _get_ear(poly, ind, eps, _i=_i+1) :
2051 // poly has no ears, look for wiskers
2052 let( wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ] )
2053 wiskers==[] ? undef : [wiskers[0]];
2054
2055
2056
2057// returns false ASA it finds some reflex vertex of poly[idxs[.]]
2058// inside the triangle different from p0 and p2
2059// note: to simplify the expressions it is assumed that the input polygon has no twists
2060function _none_inside(idxs,poly,p0,p1,p2,eps,i=0) =
2061 i>=len(idxs) ? true :
2062 let(
2063 vert = poly[idxs[i]],
2064 prev_vert = poly[select(idxs,i-1)],
2065 next_vert = poly[select(idxs,i+1)]
2066 )
2067 // check if vert prevent [p0,p1,p2] to be an ear
2068 // this conditions might have a simpler expression
2069 _tri_class([prev_vert, vert, next_vert],eps) <= 0 // reflex condition
2070 && ( // vert is a cw reflex poly vertex inside the triangle [p0,p1,p2]
2071 ( _tri_class([p0,p1,vert],eps)>0 &&
2072 _tri_class([p1,p2,vert],eps)>0 &&
2073 _tri_class([p2,p0,vert],eps)>=0 )
2074 // or it is equal to p1 and some of its adjacent edges cross the open segment (p0,p2)
2075 || ( norm(vert-p1) < eps
2076 && _is_at_left(p0,[prev_vert,p1],eps) && _is_at_left(p2,[p1,prev_vert],eps)
2077 && _is_at_left(p2,[p1,next_vert],eps) && _is_at_left(p0,[next_vert,p1],eps)
2078 )
2079 )
2080 ? false
2081 : _none_inside(idxs,poly,p0,p1,p2,eps,i=i+1);
2082
2083
2084// Function: is_polygon_clockwise()
2085// Synopsis: Determine if a 2d polygon winds clockwise.
2086// Topics: Geometry, Polygons, Clockwise
2087// See Also: clockwise_polygon(), ccw_polygon(), reverse_polygon()
2088// Usage:
2089// bool = is_polygon_clockwise(poly);
2090// Description:
2091// Return true if the given 2D simple polygon is in clockwise order, false otherwise.
2092// Results for complex (self-intersecting) polygon are indeterminate.
2093// Arguments:
2094// poly = The list of 2D path points for the perimeter of the polygon.
2095
2096// For algorithm see 2.07 here: http://www.faqs.org/faqs/graphics/algorithms-faq/
2097function is_polygon_clockwise(poly) =
2098 assert(is_path(poly,dim=2), "Input should be a 2d path")
2099 let(
2100 minx = min(poly*[1,0]),
2101 lowind = search(minx, poly, 0, 0),
2102 lowpts = select(poly,lowind),
2103 miny = min(lowpts*[0,1]),
2104 extreme_sub = search(miny, lowpts, 1, 1)[0],
2105 extreme = lowind[extreme_sub]
2106 )
2107 cross(select(poly,extreme+1)-poly[extreme],
2108 select(poly,extreme-1)-poly[extreme])<0;
2109
2110
2111// Function: clockwise_polygon()
2112// Synopsis: Return clockwise version of a polygon.
2113// Topics: Geometry, Polygons, Clockwise
2114// See Also: is_polygon_clockwise(), ccw_polygon(), reverse_polygon()
2115// Usage:
2116// newpoly = clockwise_polygon(poly);
2117// Description:
2118// Given a 2D polygon path, returns the clockwise winding version of that path.
2119// Arguments:
2120// poly = The list of 2D path points for the perimeter of the polygon.
2121function clockwise_polygon(poly) =
2122 assert(is_path(poly,dim=2), "Input should be a 2d polygon")
2123 is_polygon_clockwise(poly) ? poly : reverse_polygon(poly);
2124
2125
2126// Function: ccw_polygon()
2127// Synopsis: Return counter-clockwise version of a polygon.
2128// Topics: Geometry, Polygons, Clockwise
2129// See Also: is_polygon_clockwise(), clockwise_polygon(), reverse_polygon()
2130// Usage:
2131// newpoly = ccw_polygon(poly);
2132// Description:
2133// Given a 2D polygon poly, returns the counter-clockwise winding version of that poly.
2134// Arguments:
2135// poly = The list of 2D path points for the perimeter of the polygon.
2136function ccw_polygon(poly) =
2137 assert(is_path(poly,dim=2), "Input should be a 2d polygon")
2138 is_polygon_clockwise(poly) ? reverse_polygon(poly) : poly;
2139
2140
2141// Function: reverse_polygon()
2142// Synopsis: Reverse winding direction of polygon.
2143// Topics: Geometry, Polygons, Clockwise
2144// See Also: is_polygon_clockwise(), ccw_polygon(), clockwise_polygon()
2145// Usage:
2146// newpoly = reverse_polygon(poly)
2147// Description:
2148// Reverses a polygon's winding direction, while still using the same start point.
2149// Arguments:
2150// poly = The list of the path points for the perimeter of the polygon.
2151function reverse_polygon(poly) =
2152 let(poly=force_path(poly,"poly"))
2153 assert(is_path(poly), "Input should be a polygon")
2154 [ poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ];
2155
2156
2157// Function: reindex_polygon()
2158// Synopsis: Adjust point indexing of polygon to minimize pointwise distance to a reference polygon.
2159// Topics: Geometry, Polygons
2160// See Also: reindex_polygon(), align_polygon(), are_polygons_equal()
2161// Usage:
2162// newpoly = reindex_polygon(reference, poly);
2163// Description:
2164// Rotates and possibly reverses the point order of a 2d or 3d polygon path to optimize its pairwise point
2165// association with a reference polygon. The two polygons must have the same number of vertices and be the same dimension.
2166// The optimization is done by computing the distance, norm(reference[i]-poly[i]), between
2167// corresponding pairs of vertices of the two polygons and choosing the polygon point index rotation that
2168// makes the total sum over all pairs as small as possible. Returns the reindexed polygon. Note
2169// that the geometry of the polygon is not changed by this operation, just the labeling of its
2170// vertices. If the input polygon is 2d and is oriented opposite the reference then its point order is
2171// reversed.
2172// Arguments:
2173// reference = reference polygon path
2174// poly = input polygon to reindex
2175// Example(2D): The red dots show the 0th entry in the two input path lists. Note that the red dots are not near each other. The blue dot shows the 0th entry in the output polygon
2176// pent = subdivide_path([for(i=[0:4])[sin(72*i),cos(72*i)]],30);
2177// circ = circle($fn=30,r=2.2);
2178// reindexed = reindex_polygon(circ,pent);
2179// move_copies(concat(circ,pent)) circle(r=.1,$fn=32);
2180// color("red") move_copies([pent[0],circ[0]]) circle(r=.1,$fn=32);
2181// color("blue") translate(reindexed[0])circle(r=.1,$fn=32);
2182// Example(2D): The indexing that minimizes the total distance will not necessarily associate the nearest point of `poly` with the reference, as in this example where again the blue dot indicates the 0th entry in the reindexed result.
2183// pent = move([3.5,-1],p=subdivide_path([for(i=[0:4])[sin(72*i),cos(72*i)]],30));
2184// circ = circle($fn=30,r=2.2);
2185// reindexed = reindex_polygon(circ,pent);
2186// move_copies(concat(circ,pent)) circle(r=.1,$fn=32);
2187// color("red") move_copies([pent[0],circ[0]]) circle(r=.1,$fn=32);
2188// color("blue") translate(reindexed[0])circle(r=.1,$fn=32);
2189function reindex_polygon(reference, poly, return_error=false) =
2190 let(reference=force_path(reference,"reference"),
2191 poly=force_path(poly,"poly"))
2192 assert(is_path(reference) && is_path(poly,dim=len(reference[0])),
2193 "Invalid polygon(s) or incompatible dimensions. " )
2194 assert(len(reference)==len(poly), "The polygons must have the same length.")
2195 let(
2196 dim = len(reference[0]),
2197 N = len(reference),
2198 fixpoly = dim != 2? poly :
2199 is_polygon_clockwise(reference)
2200 ? clockwise_polygon(poly)
2201 : ccw_polygon(poly),
2202 I = [for(i=reference) 1],
2203 val = [ for(k=[0:N-1])
2204 [for(i=[0:N-1])
2205 norm(reference[i]-fixpoly[(i+k)%N]) ] ]*I,
2206 min_ind = min_index(val),
2207 optimal_poly = list_rotate(fixpoly, min_ind)
2208 )
2209 return_error? [optimal_poly, val[min_ind]] :
2210 optimal_poly;
2211
2212
2213// Function: align_polygon()
2214// Synopsis: Find best alignment of a 2d polygon to a reference 2d polygon over a set of transformations.
2215// Topics: Geometry, Polygons
2216// See Also: reindex_polygon(), align_polygon(), are_polygons_equal()
2217// Usage:
2218// newpoly = align_polygon(reference, poly, [angles], [cp], [tran], [return_ind]);
2219// Description:
2220// Find the best alignment of a specified 2D polygon with a reference 2D polygon over a set of
2221// transformations. You can specify a list or range of angles and a centerpoint or you can
2222// give a list of arbitrary 2d transformation matrices. For each transformation or angle, the polygon is
2223// reindexed, which is a costly operation so if run time is a problem, use a smaller sampling of angles or
2224// transformations. By default returns the rotated and reindexed polygon. You can also request that
2225// the best angle or the index into the transformation list be returned.
2226// Arguments:
2227// reference = reference polygon
2228// poly = polygon to rotate into alignment with the reference
2229// angles = list or range of angles to test
2230// cp = centerpoint for rotations
2231// ---
2232// tran = list of 2D transformation matrices to optimize over
2233// return_ind = if true, return the best angle (if you specified angles) or the index into tran otherwise of best alignment
2234// Example(2D): Rotating the poorly aligned light gray triangle by 105 degrees produces the best alignment, shown in blue:
2235// ellipse = yscale(3,circle(r=10, $fn=32));
2236// tri = move([-50/3,-9],
2237// subdivide_path([[0,0], [50,0], [0,27]], 32));
2238// aligned = align_polygon(ellipse,tri, [0:5:180]);
2239// color("white")stroke(tri,width=.5,closed=true);
2240// stroke(ellipse, width=.5, closed=true);
2241// color("blue")stroke(aligned,width=.5,closed=true);
2242// Example(2D,NoAxes): Translating a triangle (light gray) to the best alignment (blue)
2243// ellipse = yscale(2,circle(r=10, $fn=32));
2244// tri = subdivide_path([[0,0], [27,0], [-7,50]], 32);
2245// T = [for(x=[-10:0], y=[-30:-15]) move([x,y])];
2246// aligned = align_polygon(ellipse,tri, trans=T);
2247// color("white")stroke(tri,width=.5,closed=true);
2248// stroke(ellipse, width=.5, closed=true);
2249// color("blue")stroke(aligned,width=.5,closed=true);
2250function align_polygon(reference, poly, angles, cp, trans, return_ind=false) =
2251 let(reference=force_path(reference,"reference"),
2252 poly=force_path(poly,"poly"))
2253 assert(is_undef(trans) || (is_undef(angles) && is_undef(cp)), "Cannot give both angles/cp and trans as input")
2254 let(
2255 trans = is_def(trans) ? trans :
2256 assert( (is_vector(angles) && len(angles)>0) || valid_range(angles),
2257 "The `angle` parameter must be a range or a non void list of numbers.")
2258 [for(angle=angles) zrot(angle,cp=cp)]
2259 )
2260 assert(is_path(reference,dim=2), "reference must be a 2D polygon")
2261 assert(is_path(poly,dim=2), "poly must be a 2D polygon")
2262 assert(len(reference)==len(poly), "The polygons must have the same length.")
2263 let( // alignments is a vector of entries of the form: [polygon, error]
2264 alignments = [
2265 for(T=trans)
2266 reindex_polygon(
2267 reference,
2268 apply(T,poly),
2269 return_error=true
2270 )
2271 ],
2272 scores = column(alignments,1),
2273 minscore = min(scores),
2274 minind = [for(i=idx(scores)) if (scores[i]<minscore+EPSILON) i],
2275 dummy = is_def(angles) ? echo(best_angles = select(list(angles), minind)):0,
2276 best = minind[0]
2277 )
2278 return_ind ? (is_def(angles) ? list(angles)[best] : best)
2279 : alignments[best][0];
2280
2281
2282// Function: are_polygons_equal()
2283// Synopsis: Check if two polygons (not necessarily in the same point order) are equal.
2284// Topics: Geometry, Polygons, Comparators
2285// See Also: reindex_polygon(), align_polygon(), are_polygons_equal()
2286// Usage:
2287// bool = are_polygons_equal(poly1, poly2, [eps])
2288// Description:
2289// Returns true if poly1 and poly2 are the same polongs
2290// within given epsilon tolerance.
2291// Arguments:
2292// poly1 = first polygon
2293// poly2 = second polygon
2294// eps = tolerance for comparison
2295// Example(NORENDER):
2296// are_polygons_equal(pentagon(r=4),
2297// rot(360/5, p=pentagon(r=4))); // returns true
2298// are_polygons_equal(pentagon(r=4),
2299// rot(90, p=pentagon(r=4))); // returns false
2300function are_polygons_equal(poly1, poly2, eps=EPSILON) =
2301 let(
2302 poly1 = list_unwrap(poly1),
2303 poly2 = list_unwrap(poly2),
2304 l1 = len(poly1),
2305 l2 = len(poly2)
2306 ) l1 != l2 ? false :
2307 let( maybes = find_approx(poly1[0], poly2, eps=eps, all=true) )
2308 maybes == []? false :
2309 [for (i=maybes) if (_are_polygons_equal(poly1, poly2, eps, i)) 1] != [];
2310
2311function _are_polygons_equal(poly1, poly2, eps, st) =
2312 max([for(d=poly1-select(poly2,st,st-1)) d*d])<eps*eps;
2313
2314
2315/// Function: _is_polygon_in_list()
2316/// Topics: Polygons, Comparators
2317/// See Also: are_polygons_equal(), are_regions_equal()
2318/// Usage:
2319/// bool = _is_polygon_in_list(poly, polys);
2320/// Description:
2321/// Returns true if one of the polygons in `polys` is equivalent to the polygon `poly`.
2322/// Arguments:
2323/// poly = The polygon to search for.
2324/// polys = The list of polygons to look for the polygon in.
2325function _is_polygon_in_list(poly, polys) =
2326 ___is_polygon_in_list(poly, polys, 0);
2327
2328function ___is_polygon_in_list(poly, polys, i) =
2329 i >= len(polys)? false :
2330 are_polygons_equal(poly, polys[i])? true :
2331 ___is_polygon_in_list(poly, polys, i+1);
2332
2333
2334// Section: Convex Hull
2335
2336// This section originally based on Oskar Linde's Hull:
2337// - https://github.com/openscad/scad-utils
2338
2339
2340// Function: hull()
2341// Synopsis: Convex hull of a list of 2d or 3d points.
2342// SynTags: Ext
2343// Topics: Geometry, Hulling
2344// See Also: hull_points(), hull2d_path(), hull3d_faces()
2345// Usage:
2346// face_list_or_index_list = hull(points);
2347// Description:
2348// Takes a list of 2D or 3D points (but not both in the same list) and returns either the list of
2349// indexes into `points` that forms the 2D convex hull perimeter path, or the list of faces that
2350// form the 3d convex hull surface. Each face is a list of indexes into `points`. If the input
2351// points are co-linear, the result will be the indexes of the two extrema points. If the input
2352// points are co-planar, the results will be a simple list of vertex indices that will form a planar
2353// perimeter. Otherwise a list of faces will be returned, where each face is a simple list of
2354// vertex indices for the perimeter of the face.
2355// Arguments:
2356// points = The set of 2D or 3D points to find the hull of.
2357function hull(points) =
2358 assert(is_path(points),"Invalid input to hull")
2359 len(points[0]) == 2
2360 ? hull2d_path(points)
2361 : hull3d_faces(points);
2362
2363
2364// Module: hull_points()
2365// Synopsis: Convex hull of a list of 2d or 3d points.
2366// Topics: Geometry, Hulling
2367// See Also: hull(), hull_points(), hull2d_path(), hull3d_faces()
2368// Usage:
2369// hull_points(points, [fast]);
2370// Description:
2371// If given a list of 2D points, creates a 2D convex hull polygon that encloses all those points.
2372// If given a list of 3D points, creates a 3D polyhedron that encloses all the points. This should
2373// handle about 4000 points in slow mode. If `fast` is set to true, this should be able to handle
2374// far more. When fast mode is off, 3d hulls that lie in a plane will produce a single face of a polyhedron, which can be viewed in preview but will not render.
2375// Arguments:
2376// points = The list of points to form a hull around.
2377// fast = If true for 3d case, uses a faster cheat that may handle more points, but also may emit warnings that can stop your script if you have "Halt on first warning" enabled. Ignored for the 2d case. Default: false
2378// Example(2D):
2379// pts = [[-10,-10], [0,10], [10,10], [12,-10]];
2380// hull_points(pts);
2381// Example(3D):
2382// pts = [for (phi = [30:60:150], theta = [0:60:359]) spherical_to_xyz(10, theta, phi)];
2383// hull_points(pts);
2384module hull_points(points, fast=false) {
2385 no_children($children);
2386 check = assert(is_path(points))
2387 assert(len(points)>=3, "Point list must contain 3 points");
2388 attachable(){
2389 if (len(points[0])==2)
2390 hull() polygon(points=points);
2391 else if (len(points)==3)
2392 polyhedron(points=points, faces=[[0,1,2]]);
2393 else {
2394 if (fast) {
2395 extra = len(points)%3;
2396 faces = [
2397 [for(i=[0:1:extra+2])i], // If vertex count not divisible by 3, combine extras with first 3
2398 for(i=[extra+3:3:len(points)-3])[i,i+1,i+2]
2399 ];
2400 hull() polyhedron(points=points, faces=faces);
2401 } else {
2402 faces = hull(points);
2403 if (is_num(faces[0])){
2404 if (len(faces)<=2) echo("Hull contains only two points");
2405 else polyhedron(points=points, faces=[faces]);
2406 }
2407 else polyhedron(points=points, faces=faces);
2408 }
2409 }
2410 union();
2411 }
2412}
2413
2414
2415
2416function _backtracking(i,points,h,t,m,all) =
2417 m<t || _is_cw(points[i], points[h[m-1]], points[h[m-2]],all) ? m :
2418 _backtracking(i,points,h,t,m-1,all) ;
2419
2420// clockwise check (2d)
2421function _is_cw(a,b,c,all) =
2422 all ? cross(a-c,b-c)<=EPSILON*norm(a-c)*norm(b-c) :
2423 cross(a-c,b-c)<-EPSILON*norm(a-c)*norm(b-c);
2424
2425
2426// Function: hull2d_path()
2427// Synopsis: Convex hull of a list of 2d points.
2428// Topics: Geometry, Hulling
2429// See Also: hull(), hull_points(), hull2d_path(), hull3d_faces()
2430// Usage:
2431// index_list = hull2d_path(points,all)
2432// Description:
2433// Takes a list of arbitrary 2D points, and finds the convex hull polygon to enclose them.
2434// Returns a path as a list of indices into `points`.
2435// When all==true, returns extra points that are on edges of the hull.
2436// Arguments:
2437// points = list of 2d points to get the hull of.
2438// all = when true, includes all points on the edges of the convex hull. Default: false.
2439// Example(2D):
2440// pts = [[-10,-10], [0,10], [10,10], [12,-10]];
2441// path = hull2d_path(pts);
2442// move_copies(pts) color("red") circle(1,$fn=12);
2443// polygon(points=pts, paths=[path]);
2444//
2445// Code based on this method:
2446// https://www.hackerearth.com/practice/math/geometry/line-sweep-technique/tutorial/
2447//
2448function hull2d_path(points, all=false) =
2449 assert(is_path(points,2),"Invalid input to hull2d_path")
2450 len(points) < 2 ? [] :
2451 let( n = len(points),
2452 ip = sortidx(points) )
2453 // lower hull points
2454 let( lh =
2455 [ for( i = 2,
2456 k = 2,
2457 h = [ip[0],ip[1]]; // current list of hull point indices
2458 i <= n;
2459 k = i<n ? _backtracking(ip[i],points,h,2,k,all)+1 : k,
2460 h = i<n ? [for(j=[0:1:k-2]) h[j], ip[i]] : [],
2461 i = i+1
2462 ) if( i==n ) h ][0] )
2463 // concat lower hull points with upper hull ones
2464 [ for( i = n-2,
2465 k = len(lh),
2466 t = k+1,
2467 h = lh; // current list of hull point indices
2468 i >= -1;
2469 k = i>=0 ? _backtracking(ip[i],points,h,t,k,all)+1 : k,
2470 h = [for(j=[0:1:k-2]) h[j], if(i>0) ip[i]],
2471 i = i-1
2472 ) if( i==-1 ) h ][0] ;
2473
2474
2475function _hull_collinear(points) =
2476 let(
2477 a = points[0],
2478 i = max_index([for(pt=points) norm(pt-a)]),
2479 n = points[i] - a
2480 )
2481 norm(n)==0 ? [0]
2482 :
2483 let(
2484 points1d = [ for(p = points) (p-a)*n ],
2485 min_i = min_index(points1d),
2486 max_i = max_index(points1d)
2487 ) [min_i, max_i];
2488
2489
2490
2491// Function: hull3d_faces()
2492// Synopsis: Convex hull of a list of 3d points.
2493// Topics: Geometry, Hulling
2494// See Also: hull(), hull_points(), hull2d_path(), hull3d_faces()
2495// Usage:
2496// faces = hull3d_faces(points)
2497// Description:
2498// Takes a list of arbitrary 3D points, and finds the convex hull polyhedron to enclose
2499// them. Returns a list of triangular faces, where each face is a list of indexes into the given `points`
2500// list. The output will be valid for use with the polyhedron command, but may include vertices that are in the interior of a face of the hull, so it is not
2501// necessarily the minimal representation of the hull.
2502// If all points passed to it are coplanar, then the return is the list of indices of points
2503// forming the convex hull polygon.
2504// Example(3D):
2505// pts = [[-20,-20,0], [20,-20,0], [0,20,5], [0,0,20]];
2506// faces = hull3d_faces(pts);
2507// move_copies(pts) color("red") sphere(1);
2508// %polyhedron(points=pts, faces=faces);
2509function hull3d_faces(points) =
2510 assert(is_path(points,3),"Invalid input to hull3d_faces")
2511 len(points) < 3 ? count(len(points))
2512 : let ( // start with a single non-collinear triangle
2513 tri = _noncollinear_triple(points, error=false)
2514 )
2515 tri==[] ? _hull_collinear(points)
2516 : let(
2517 a = tri[0],
2518 b = tri[1],
2519 c = tri[2],
2520 plane = plane3pt_indexed(points, a, b, c),
2521 d = _find_first_noncoplanar(plane, points)
2522 )
2523 d == len(points)
2524 ? /* all coplanar*/
2525 let (
2526 pts2d = project_plane([points[a], points[b], points[c]],points),
2527 hull2d = hull2d_path(pts2d)
2528 ) hull2d
2529 : let(
2530 remaining = [for (i = [0:1:len(points)-1]) if (i!=a && i!=b && i!=c && i!=d) i],
2531 // Build an initial tetrahedron.
2532 // Swap b, c if d is in front of triangle t.
2533 ifop = _is_point_above_plane(plane, points[d]),
2534 bc = ifop? [c,b] : [b,c],
2535 b = bc[0],
2536 c = bc[1],
2537 triangles = [
2538 [a,b,c],
2539 [d,b,a],
2540 [c,d,a],
2541 [b,d,c]
2542 ],
2543 // calculate the plane equations
2544 planes = [ for (t = triangles) plane3pt_indexed(points, t[0], t[1], t[2]) ]
2545 ) _hull3d_iterative(points, triangles, planes, remaining);
2546
2547
2548// Adds the remaining points one by one to the convex hull
2549function _hull3d_iterative(points, triangles, planes, remaining, _i=0) = //let( EPSILON=1e-12 )
2550 _i >= len(remaining) ? triangles :
2551 let (
2552 // pick a point
2553 i = remaining[_i],
2554 // evaluate the triangle plane equations at point i
2555 planeq_val = planes*[each points[i], -1],
2556 // find the triangles that are in conflict with the point (point not inside)
2557 conflicts = [for (i = [0:1:len(planeq_val)-1]) if (planeq_val[i]>EPSILON) i ],
2558 // collect the halfedges of all triangles that are in conflict
2559 halfedges = [
2560 for(c = conflicts, i = [0:2])
2561 [triangles[c][i], triangles[c][(i+1)%3]]
2562 ],
2563 // find the outer perimeter of the set of conflicting triangles
2564 horizon = _remove_internal_edges(halfedges),
2565 // generate new triangles connecting point i to each horizon halfedge vertices
2566 tri2add = [ for (h = horizon) concat(h,i) ],
2567 // add tria2add and remove conflict triangles
2568 new_triangles =
2569 concat( tri2add,
2570 [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) triangles[i] ]
2571 ),
2572 // add the plane equations of new added triangles and remove the plane equations of the conflict ones
2573 new_planes =
2574 [ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) ,
2575 for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ]
2576 ) _hull3d_iterative(
2577 points,
2578 new_triangles,
2579 new_planes,
2580 remaining,
2581 _i+1
2582 );
2583
2584
2585function _remove_internal_edges(halfedges) = [
2586 for (h = halfedges)
2587 if (!in_list(reverse(h), halfedges))
2588 h
2589];
2590
2591function _find_first_noncoplanar(plane, points, i=0) =
2592 (i >= len(points) || !are_points_on_plane([points[i]],plane))? i :
2593 _find_first_noncoplanar(plane, points, i+1);
2594
2595
2596// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
2597
2598
2599
2600
2601
2602// Section: Convex Sets
2603
2604
2605// Function: is_polygon_convex()
2606// Synopsis: Check if a polygon is convex.
2607// Topics: Geometry, Convexity, Test
2608// See Also: clockwise_polygon(), ccw_polygon()
2609// Usage:
2610// bool = is_polygon_convex(poly, [eps]);
2611// Description:
2612// Returns true if the given 2D or 3D polygon is convex.
2613// The result is meaningless if the polygon is not simple (self-crossing) or non coplanar.
2614// If the points are collinear or not coplanar an error may be generated.
2615// Arguments:
2616// poly = Polygon to check.
2617// eps = Tolerance for the collinearity and coplanarity tests. Default: EPSILON.
2618// Example:
2619// test1 = is_polygon_convex(circle(d=50)); // Returns: true
2620// test2 = is_polygon_convex(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
2621// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
2622// test = is_polygon_convex(spiral); // Returns: false
2623function is_polygon_convex(poly,eps=EPSILON) =
2624 assert(is_path(poly), "The input should be a 2D or 3D polygon." )
2625 let(
2626 lp = len(poly),
2627 p0 = poly[0]
2628 )
2629 assert( lp>=3 , "A polygon must have at least 3 points" )
2630 let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] )
2631 len(p0)==2
2632 ? let( size = max([for(p=poly) norm(p-p0)]), tol=pow(size,2)*eps )
2633 assert( size>eps, "The polygon is self-crossing or its points are collinear" )
2634 min(crosses) >=-tol || max(crosses)<=tol
2635 : let( ip = _noncollinear_triple(poly,error=false,eps=eps) )
2636 assert( ip!=[], "The points are collinear")
2637 let(
2638 crx = cross(poly[ip[1]]-poly[ip[0]],poly[ip[2]]-poly[ip[1]]),
2639 nrm = crx/norm(crx),
2640 plane = concat(nrm, nrm*poly[0]),
2641 prod = crosses*nrm,
2642 size = norm(poly[ip[1]]-poly[ip[0]]),
2643 tol = pow(size,2)*eps
2644 )
2645 assert(_pointlist_greatest_distance(poly,plane) < size*eps, "The polygon points are not coplanar")
2646 let(
2647 minc = min(prod),
2648 maxc = max(prod) )
2649 minc>=-tol || maxc<=tol;
2650
2651
2652// Function: convex_distance()
2653// Synopsis: Compute distance between convex hull of two point lists.
2654// Topics: Geometry, Convexity, Distance
2655// See also: convex_collision(), hull()
2656// Usage:
2657// dist = convex_distance(points1, points2,eps);
2658// Description:
2659// Returns the smallest distance between a point in convex hull of `points1`
2660// and a point in the convex hull of `points2`. All the points in the lists
2661// should have the same dimension, either 2D or 3D.
2662// A zero result means the hulls intercept whithin a tolerance `eps`.
2663// Arguments:
2664// points1 = first list of 2d or 3d points.
2665// points2 = second list of 2d or 3d points.
2666// eps = tolerance in distance evaluations. Default: EPSILON.
2667// Example(2D):
2668// pts1 = move([-3,0], p=square(3,center=true));
2669// pts2 = rot(a=45, p=square(2,center=true));
2670// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
2671// polygon(pts1);
2672// polygon(pts2);
2673// polygon(pts3);
2674// echo(convex_distance(pts1,pts2)); // Returns: 0.0857864
2675// echo(convex_distance(pts2,pts3)); // Returns: 0
2676// Example(3D):
2677// sphr1 = sphere(2,$fn=10);
2678// sphr2 = move([4,0,0], p=sphr1);
2679// sphr3 = move([4.5,0,0], p=sphr1);
2680// vnf_polyhedron(sphr1);
2681// vnf_polyhedron(sphr2);
2682// echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0
2683// echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5
2684function convex_distance(points1, points2, eps=EPSILON) =
2685 assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])),
2686 "The input lists should be compatible consistent non empty lists of points.")
2687 assert(len(points1[0])==2 || len(points1[0])==3 ,
2688 "The input points should be 2d or 3d points.")
2689 let( d = points1[0]-points2[0] )
2690 norm(d)<eps ? 0 :
2691 let( v = _support_diff(points1,points2,-d) )
2692 norm(_GJK_distance(points1, points2, eps, 0, v, [v]));
2693
2694
2695// Finds the vector difference between the hulls of the two pointsets by the GJK algorithm
2696// Based on:
2697// http://www.dtecta.com/papers/jgt98convex.pdf
2698function _GJK_distance(points1, points2, eps=EPSILON, lbd, d, simplex=[]) =
2699 let( nrd = norm(d) ) // distance upper bound
2700 nrd<eps ? d :
2701 let(
2702 v = _support_diff(points1,points2,-d),
2703 lbd = max(lbd, d*v/nrd), // distance lower bound
2704 close = (nrd-lbd <= eps*nrd)
2705 )
2706 close ? d :
2707 let( newsplx = _closest_simplex(concat(simplex,[v]),eps) )
2708 _GJK_distance(points1, points2, eps, lbd, newsplx[0], newsplx[1]);
2709
2710
2711// Function: convex_collision()
2712// Synopsis: Check whether the convex hulls of two point lists intersect.
2713// Topics: Geometry, Convexity, Collision, Intersection
2714// See also:
2715// convex_distance(), hull()
2716// Usage:
2717// bool = convex_collision(points1, points2, [eps]);
2718// Description:
2719// Returns `true` if the convex hull of `points1` intersects the convex hull of `points2`
2720// otherwise, `false`.
2721// All the points in the lists should have the same dimension, either 2D or 3D.
2722// This function is tipically faster than `convex_distance` to find a non-collision.
2723// Arguments:
2724// points1 = first list of 2d or 3d points.
2725// points2 = second list of 2d or 3d points.
2726// eps - tolerance for the intersection tests. Default: EPSILON.
2727// Example(2D):
2728// pts1 = move([-3,0], p=square(3,center=true));
2729// pts2 = rot(a=45, p=square(2,center=true));
2730// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
2731// polygon(pts1);
2732// polygon(pts2);
2733// polygon(pts3);
2734// echo(convex_collision(pts1,pts2)); // Returns: false
2735// echo(convex_collision(pts2,pts3)); // Returns: true
2736// Example(3D):
2737// sphr1 = sphere(2,$fn=10);
2738// sphr2 = move([4,0,0], p=sphr1);
2739// sphr3 = move([4.5,0,0], p=sphr1);
2740// vnf_polyhedron(sphr1);
2741// vnf_polyhedron(sphr2);
2742// echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true
2743// echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false
2744//
2745function convex_collision(points1, points2, eps=EPSILON) =
2746 assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])),
2747 "The input lists should be compatible consistent non empty lists of points.")
2748 assert(len(points1[0])==2 || len(points1[0])==3 ,
2749 "The input points should be 2d or 3d points.")
2750 let( d = points1[0]-points2[0] )
2751 norm(d)<eps ? true :
2752 let( v = _support_diff(points1,points2,-d) )
2753 _GJK_collide(points1, points2, v, [v], eps);
2754
2755
2756// Based on the GJK collision algorithms found in:
2757// http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
2758// or
2759// http://www.dtecta.com/papers/jgt98convex.pdf
2760function _GJK_collide(points1, points2, d, simplex, eps=EPSILON) =
2761 norm(d) < eps ? true : // does collide
2762 let( v = _support_diff(points1,points2,-d) )
2763 v*d > eps*eps ? false : // no collision
2764 let( newsplx = _closest_simplex(concat(simplex,[v]),eps) )
2765 norm(v-newsplx[0])<eps ? norm(v)<eps :
2766 _GJK_collide(points1, points2, newsplx[0], newsplx[1], eps);
2767
2768
2769// given a simplex s, returns a pair:
2770// - the point of the s closest to the origin
2771// - the smallest sub-simplex of s that contains that point
2772function _closest_simplex(s,eps=EPSILON) =
2773 len(s)==2 ? _closest_s1(s,eps) :
2774 len(s)==3 ? _closest_s2(s,eps) :
2775 len(s)==4 ? _closest_s3(s,eps) :
2776 assert(false, "Internal error.");
2777
2778
2779// find the point of a 1-simplex closest to the origin
2780function _closest_s1(s,eps=EPSILON) =
2781 norm(s[1]-s[0])<=eps*(norm(s[0])+norm(s[1]))/2 ? [ s[0], [s[0]] ] :
2782 let(
2783 c = s[1]-s[0],
2784 t = -s[0]*c/(c*c)
2785 )
2786 t<0 ? [ s[0], [s[0]] ] :
2787 t>1 ? [ s[1], [s[1]] ] :
2788 [ s[0]+t*c, s ];
2789
2790
2791// find the point of a 2-simplex closest to the origin
2792function _closest_s2(s, eps=EPSILON) =
2793 // considering that s[2] was the last inserted vertex in s by GJK,
2794 // the plane orthogonal to the triangle [ origin, s[0], s[1] ] that
2795 // contains [s[0],s[1]] have the origin and s[2] on the same side;
2796 // that reduces the cases to test and the only possible simplex
2797 // outcomes are s, [s[0],s[2]] and [s[1],s[2]]
2798 let(
2799 area = cross(s[2]-s[0], s[1]-s[0]),
2800 area2 = area*area // tri area squared
2801 )
2802 area2<=eps*max([for(si=s) pow(si*si,2)]) // degenerate tri
2803 ? norm(s[2]-s[0]) < norm(s[2]-s[1])
2804 ? _closest_s1([s[1],s[2]])
2805 : _closest_s1([s[0],s[2]])
2806 : let(
2807 crx1 = cross(s[0], s[2])*area,
2808 crx2 = cross(s[1], s[0])*area,
2809 crx0 = cross(s[2], s[1])*area
2810 )
2811 // all have the same signal -> origin projects inside the tri
2812 max(crx1, crx0, crx2) < 0 || min(crx1, crx0, crx2) > 0
2813 ? // baricentric coords of projection
2814 [ [abs(crx0),abs(crx1),abs(crx2)]*s/area2, s ]
2815 : let(
2816 cl12 = _closest_s1([s[1],s[2]]),
2817 cl02 = _closest_s1([s[0],s[2]])
2818 )
2819 norm(cl12[0])<norm(cl02[0]) ? cl12 : cl02;
2820
2821
2822// find the point of a 3-simplex closest to the origin
2823function _closest_s3(s,eps=EPSILON) =
2824 let( nr = cross(s[1]-s[0],s[2]-s[0]),
2825 sz = [ norm(s[0]-s[1]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] )
2826 norm(nr)<=eps*pow(max(sz),2)
2827 ? let( i = max_index(sz) )
2828 _closest_s2([ s[i], s[(i+1)%3], s[3] ], eps) // degenerate case
2829 : // considering that s[3] was the last inserted vertex in s by GJK,
2830 // the only possible outcomes will be:
2831 // s or some of the 3 faces of s containing s[3]
2832 let(
2833 tris = [ [s[0], s[1], s[3]],
2834 [s[1], s[2], s[3]],
2835 [s[2], s[0], s[3]] ],
2836 cntr = sum(s)/4,
2837 // indicator of the tris facing the origin
2838 facing = [for(i=[0:2])
2839 let( nrm = _tri_normal(tris[i]) )
2840 if( ((nrm*(s[i]-cntr))>0)==(nrm*s[i]<0) ) i ]
2841 )
2842 len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex
2843 len(facing)==1 ? _closest_s2(tris[facing[0]], eps) :
2844 let( // look for the origin-facing tri closest to the origin
2845 closest = [for(i=facing) _closest_s2(tris[i], eps) ],
2846 dist = [for(cl=closest) norm(cl[0]) ],
2847 nearest = min_index(dist)
2848 )
2849 closest[nearest];
2850
2851
2852function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]);
2853
2854
2855function _support_diff(p1,p2,d) =
2856 let( p1d = p1*d, p2d = p2*d )
2857 p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]];
2858
2859
2860// Section: Rotation Decoding
2861
2862// Function: rot_decode()
2863// Synopsis: Extract axis and rotation angle from a rotation matrix.
2864// Topics: Affine, Matrices, Transforms
2865// Usage:
2866// info = rot_decode(rotation,[long]); // Returns: [angle,axis,cp,translation]
2867// Description:
2868// Given an input 3D rigid transformation operator (one composed of just rotations and translations) represented
2869// as a 4x4 matrix, compute the rotation and translation parameters of the operator. Returns a list of the
2870// four parameters, the angle, in the interval [0,180], the rotation axis as a unit vector, a centerpoint for
2871// the rotation, and a translation. If you set `parms = rot_decode(rotation)` then the transformation can be
2872// reconstructed from parms as `move(parms[3]) * rot(a=parms[0],v=parms[1],cp=parms[2])`. This decomposition
2873// makes it possible to perform interpolation. If you construct a transformation using `rot` the decoding
2874// may flip the axis (if you gave an angle outside of [0,180]). The returned axis will be a unit vector, and
2875// the centerpoint lies on the plane through the origin that is perpendicular to the axis. It may be different
2876// than the centerpoint you used to construct the transformation.
2877// .
2878// If you set `long` to true then return the reversed rotation, with the angle in [180,360].
2879// Arguments:
2880// rotation = rigid transformation to decode
2881// long = if true return the "long way" around, with the angle in [180,360]. Default: false
2882// Example:
2883// info = rot_decode(rot(45));
2884// // Returns: [45, [0,0,1], [0,0,0], [0,0,0]]
2885// info = rot_decode(rot(a=37, v=[1,2,3], cp=[4,3,-7])));
2886// // Returns: [37, [0.26, 0.53, 0.80], [4.8, 4.6, -4.6], [0,0,0]]
2887// info = rot_decode(left(12)*xrot(-33));
2888// // Returns: [33, [-1,0,0], [0,0,0], [-12,0,0]]
2889// info = rot_decode(translate([3,4,5]));
2890// // Returns: [0, [0,0,1], [0,0,0], [3,4,5]]
2891function rot_decode(M,long=false) =
2892 assert(is_matrix(M,4,4) && approx(M[3],[0,0,0,1]), "Input matrix must be a 4x4 matrix representing a 3d transformation")
2893 let(R = submatrix(M,[0:2],[0:2]))
2894 assert(approx(det3(R),1) && approx(norm_fro(R * transpose(R)-ident(3)),0),"Input matrix is not a rotation")
2895 let(
2896 translation = [for(row=[0:2]) M[row][3]], // translation vector
2897 largest = max_index([R[0][0], R[1][1], R[2][2]]),
2898 axis_matrix = R + transpose(R) - (matrix_trace(R)-1)*ident(3), // Each row is on the rotational axis
2899 // Construct quaternion q = c * [x sin(theta/2), y sin(theta/2), z sin(theta/2), cos(theta/2)]
2900 q_im = axis_matrix[largest],
2901 q_re = R[(largest+2)%3][(largest+1)%3] - R[(largest+1)%3][(largest+2)%3],
2902 c_sin = norm(q_im), // c * sin(theta/2) for some c
2903 c_cos = abs(q_re) // c * cos(theta/2)
2904 )
2905 approx(c_sin,0) ? [0,[0,0,1],[0,0,0],translation] :
2906 let(
2907 angle = 2*atan2(c_sin, c_cos), // This is supposed to be more accurate than acos or asin
2908 axis = (q_re>=0 ? 1:-1)*q_im/c_sin,
2909 tproj = translation - (translation*axis)*axis, // Translation perpendicular to axis determines centerpoint
2910 cp = (tproj + cross(axis,tproj)*c_cos/c_sin)/2
2911 )
2912 [long ? 360-angle:angle,
2913 long? -axis : axis,
2914 cp,
2915 (translation*axis)*axis];
2916
2917
2918
2919
2920// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap